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I. INTRODUCTION 

The Canonical Function Method (CFM) is a power- 
ful means for solving the Radial Schrodinger Equation 
(RSE), a singular boundary value problem (SBVP) with 
a stringent dual requirement of regularity near the ori- 
gin (r ~ 0) where the potential is large and near infinity 
(r — > oo) where the potential is very small. 

The CFM turns this SBVP into a regular initial value 
problem and allows full and accurate determination of 
the spectrum of the Schrodinger operator. 

The CFM can handle bound states and scattering 
states, spinless and spin-dependent, single channel or 
multi-channel problems and possesses several features 



that makes it of interest for solving a large variety of 
problems in Quantum Mechanics of Atoms, Molecules 
and Scattering. Those capabilities are based essentially 
on the following main characteristics: 

• Its capability to transform a SVBP into a regular 
initial value problem (PJVP), means mathemati- 
cally that the CFM belongs to a family of Invari- 
ant Embedding Methods (IEM). Those essentially 
transform a boundary value problem into an initial 
value problem. The CFM is even more powerful 
than standard IEM techniques since it turns the 
nature of the problem from singular to regular. 

• The evaluation of the Schrodinger operator spec- 
trum is done without performing diagonalization, 
bypassing the evaluation of the eigenfunctions. 
This allows to preserve a high degree of numerical 
precision that is required in solving sensitive eigen- 
value problems arising in Cold Molecules or large 
n-limit problems where levels to evaluate are very 
close to the continuum limit or dissociation. 

The numerical precision gained with the bypass of in- 
termediate digonalization operations is reminiscent of the 
Golub-Reinsch algorithm (see for instance ref. (f50h ) used 
for the singular value decomposition of arbitrary rectan- 
gular matrices. 

This review is organised as follows: The next section 
is a description of the CFM with highlights of its mathe- 
matical aspects and its generalisation to the multichannel 
case. In section III the mathematical and numerical tech- 
niques of the radial Schrodinger equation are discussed 
and the full procedure for dealing with general potentials, 
energy levels and evaluation of phaseshifts is detailed in 
section IV. In the rest of the sections, we apply the CFM 
method systematically to a variety of problems highlight- 
ing its successes and showing its limitations whenever 
they show up. While, the versatility of the method is 
illustrated with the help of this selection of examples, 
the list below is by no means exhaustive but serves the 
purpose of underlining the breadth of applications of the 
CFM. 
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• Energy levels for regular and singular potentials 
(radial case) 

• Potential estimation (in parametric or pseudo 
form) from spectroscopic data 

• Test of the Bohr Correspondence Principle 

• Vibrational energy levels of Cold Molecules and ap- 
plication to the 23 Na2 molecule in the 0~ and l u 
electronic states. The Lennard- Jones molecule case 
is also studied. 

• Local and Non-Local Exchange problems. 

• Accurate Phase shift evaluation for regular and sin- 
gular potentials. 

In the last section before the conclusion, we treat the 
general but non obvious (regular and singular) ID poten- 
tials to illustrate how the method can be applied in this 
case of lesser interest than the 3D's. 

In the Appendix we provide information linking the 
various units used since the review contains results span- 
ning several fields of physics where different systems of 
units are used. The boundary conditions affecting the 
spectrum determination via the matching conditions at 
a single radial point are also discussed in the Appendix. 



II. THE CANONICAL FUNCTION METHOD (CFM) 

The CFM (|29l ) is a powerful means for solving the 
Radial Schrodinger Equation (RSE). The mathematical 
difficulty of the RSE lies in the fact it is a singular 
boundary value problem. The CFM turns it into a 
regular initial value problem and allows the full deter- 
mination of the spectrum of the Schrodinger operator 
bypassing the evaluation of the eigenfunctions. 

The partial wave form of the RSE is written as: 



Hf d 2 ui{E-r) 
2fi dr 2 



h 2 i(i + 1) 



lim ui(r) = 0; lim ui(r) = 



(2) 



The CFM developed initially by Kobeissi (|35T ) and his 
coworkers to integrate the RSE, consists of writing the 
general solution y(r) representing the probability ampli- 
tude ui(E;r) as a function of the radial distance r in 
terms of two basis functions a(E; r) and f3(E; r) for some 
energy E. 



Generally, the RSE is rewritten in a system of units such 
that h = 1, 2fi — 1 (see Appendix on units): 



d 2 y(r) 
dr 2 



V(r) 



1(1 + 1) 



- E 



y(r) 



(3) 



At a selected distance r , a well defined set of initial 
conditions are satisfied by the canonical functions and 
their derivatives ie: a(E;ro) = 1 with a'(E;ro) = and 
P(E;r ) = with &'{E;r ) = 1. Thus we write: 



y(r) = y(r )a(E; r) + y / (r )P(E; r) 



(4) 



The method of solving the RSE is to proceed from r$ 
simultaneously towards the origin (r — > 0) and towards 
infinity (r — ► oo). Expressing the continuity condition of 
the " wavefunction" y(r) and its derivative y'(r) at the 
point ro from the left (r — > 0) and right (r — > oo) yields 
(see Appendix on matching and boundary conditions): 



y'(r a ) 
y(ro) 
y'(ro) 



y(ro) 



_a(E-0) 

a(E; oo) 
'(3(E;oo) 



(•5) 



When the integration is performed, the ratio of the r 
dependent canonical functions is monitored until satura- 
tion with respect to r is reached at both limits (r — > 
and r — > oo). The saturation of the "f"!!' r ? ratio with r 



(3(E-r) 



yields a position independent eigenvalue function F(E). 
The latter is mathematically defined with the help of two 
associated energy functions: 



I- (E) = lim 



a(E;r) 
"/?(£; r) 



+ + — „ 2 ME;r) = Eui(E;r) 

(1) 

where \i is the reduced mass and Ui(E; r) is the reduced 
probability amplitude for orbital angular momentum I 
and eigenvalue E. 

The boundary conditions are: 



and: 



as: 



l + (E) 



lim 



a(E;r) 



F(E) = l + (E)-l_(E) = 



~y'{ro)~ 




|Y(ro)l 


. yi r o) _ 


+ 


. y( r o) . 



(0) 



(7) 



(8) 



its zeroes expressing the continuity of y(r) and its 
derivative y'(r) at the point ro (matching conditions) 
yield the spectrum of the RSE. In practice we have: 



and: 



I -(E) 



l + (E) 



a(E; r n 



0(E;r min ) 



a(E;r max ) 
(3(E; r max ) 



(9) 



(10) 
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with r m i n , f max the radial coordinates where satu- 
ration is observed (within some predefined tolerance) 
respectively in 1-{E) and l+(E). An example of typical 
behaviour of F(E) is displayed in Fig. [T] The eigen- 
functions may be obtained for any E — Ej~ where E% is 
a zero of F(E). 



pression of the total potential is defined as: 
1(1 + 1 



V(r) 



= Y,ln P) X n ,X£lp (11) 
rei P n= Q 



In order to perform integration, we use the following 
variable step difference equation: 



y P+ i=y p + h p y' p + J2ck p) K (12) 

n=2 




-4 



100 200 300 400 500 600 700 800 900 1000 
E 



FIG. 1 Typical behavior of the eigenvalue function with en- 
ergy. The vertical lines indicate the eigenvalue position. The 
graph of the eigenvalue function has an approximate tan(_B) 
shape versus the energy E. 

The eigenvalue function definition depends on the 
type of boundary conditions at hand. In the Appendix 
we describe the general boundary conditions case and 
the corresponding matching conditions along with the 
corresponding eigenvalue function. 

The t&n(E) shape of F(E) provides a deep insight into 
the physical significance of the CFM method. The latter 
transforms a SBVP from the open interval [0, oof to the 
finite interval [r m i n ,r max ] defined by the saturation co- 
ordinates of the energy functions. This means the CFM 
maps an arbitrary potential V(r) onto the infinite square 
well problem defined by V(r) = in the open interval 
]r m in,r m ax[ and V(r min ) = V(r max ) = oo (see section 
VII for additional information). In fact, for the infi- 
nite square well problem F(E) has a t&n(E) (see also 
ref. (HI)). 

The speed and accuracy of the CFM method have been 
tested and compared to standard integration algorithms 
such as the order four Runge-Kutta (RK4) method, Nu- 
merov etc.. in a variety of cases and for a wide of range 
of potentials. 

More specifically, tq being selected as the starting point 
for the integration, the r-axis is divided into intervals 
Ip = [r pi r p+ i] in a way such that t he p otential can be 
series expanded over I p (see refs (fill fl7n). 

(v) 

The expansion uses coefficients 7„ so that the local ex- 



y' P+1 =y' P + J2 nC n )h r 1 ( 13 ) 

n=2 

where the y p ~ y(r) are a particular set of canoni- 
cal functions and the y' p — y'(r) are their derivatives in 

the interval [r p ,r p+ i]. The C n p ^ are given by the recur- 
sion formula for energy E obtained from series expanding 
eq.[3] 

n 

(n + 2)(n + l)C n % = C^l^L ~ EC™ (14) 

with C (p) = y p and C {p) = y' p . 

In practice, the above sums over n are truncated to 
a cutoff value N (instead of oo) chosen large enough 
so that the remainder error is less than a selected 
truncation error et ■ 

Given N and er, the integration step h p for the interval 

j_ 

[r p ,r p +i] is deduced from h p — (t\,/\C n P \) N ■ 
At the starting point of integration ro, the starting step 
size ho is thus determined, and the procedure is repeated 
at n = ro + ho leading to hi and so forth until we reach 
r P +i = r p + h p after p steps. 

These equations allow the propagation of the solution 
from one point to the next using a variable step and 
a local series expansion of variable order controlled by 
the truncation error ex- The accuracy of the results are 
monitored with respect to the decrease of e$ and the 
procedure is referred to as the VSCA (Variable Step with 
Controlled Accuracy) method. 

Generally, one avoids using or evaluating the radial 
wavefunction but if one insists on evaluating it, the 
canonical functions a(E;r) and f3(E;r) are used to de- 
termine the radial wavefunction at any energy via the 
expression: 

ui{E-r) = u l (E;r )a(E;r)+u' l (E;r )P(E;r) (15) 

where ui(E;r) and u[(E;ro) are the radial wavefunction 
and its derivative at the initial distance ro- 

Before proceeding any further, the essential test of any 
method that pretends solving the RSE is the Coulomb 
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potential and the essential case to test the accuracy and 
reliability of the CFM is the (textbook) Hydrogen atom. 

The CFM results are shown in Table. U along with the 
exact analytical results. 



Index CFM (Ry) Exact (Ry) 

1 -1.00000 -1.00000 

2 -0.250000 -0.250000 



3 


-0.111111 


-0.111111 


4 


-6.25000( 


-2) 


-6.25000( 


-2) 


5 


-4.00000( 


-2) 


-4.00000( 


-2) 


6 


-2.77778( 


-2) 


-2.77778( 


-2) 


7 


-2.04082( 


-2) 


-2.04082( 


-2) 


8 


-1.56250( 


-2) 


-1.56250( 


-2) 


9 


-1.23457( 


-2) 


-1.23457( 


-2) 


10 


-1.000001 


-2) 


-1.00000( 


-2) 


11 


-8.26446( 


-3) 


-8.26446( 


-3) 


12 


-6.94444( 


-3) 


-6.94444( 


-3) 


13 


-5.91716( 


-3) 


-5.91716( 


-3) 


14 


-5.10204( 


-3) 


-5.10204( 


-3) 


15 


-4.44445( 


-3) 


-4.44445( 


-3) 


16 


-3.90625( 


-3) 


-3.90625( 


-3) 


17 


-3.46021( 


-3) 


-3.46021( 


-3) 


18 


-3.08642( 


-3) 


-3.08642( 


-3) 


19 


-2.77008( 


-3) 


-2.77008( 


-3) 


20 


-2.50000( 


-3) 


-2.50000( 


-3) 


21 


-2.26757( 


-3) 


-2.26757( 


-3) 


22 


-2.06612( 


-3) 


-2.06612( 


-3) 


23 


-1.89036( 


-3) 


-1.89036( 


-3) 


24 


-1.7361K 


-3) 


-1.7361K 


-3) 



TABLE I Energy levels of the Hydrogen atom. Middle col- 
umn values are the CFM results whereas the last column val- 
ues are the corresponding exact analytically obtained values. 
The numbers in parenthesis represent a power of 10. 

It is remarkable to notice that all digits (calculated by 
CFM and analytically) are all same. 

The classical Morse potential is the simplest model 
for the evaluation of vibrational spectra of diatomic 
molecules. First of all, we have analytical expressions 
for the levels, secondly it provides a stringent test of 
the CFM before embarking into more sophisticated cases 
such as weakly bound or cold molecules (also called long- 
range molecules). 

The Morse potential is given by: 

V(r) = D[l — exp(— a{r — r e })] 2 - D (16) 

with the values D, a, r e equal respectively to 188.4355, 
0.711248, 1.9975 in au. The analytic expression for the 
levels is: 



a 2f^ 

£„ = - — (v/p/a-n-1/2) 2 , (17) 

with max n < ^/2fiD/a — 1/2. Hence the number of 
levels is given by: N — y/2jlD/a — 1/2. 

Working with units such that Ti = 1 and 2/i = 1, the 
Morse potential and the eigenvalue function F(E) are 
displayed in the figures below. The table contains the 
levels calculated by CFM and compared to the analyti- 
cal case. Again, like in the pure Coulomb case treated 
previously, the agreement is perfect and we found all the 
levels (N = 19) as predicted analytically. 

Morse potential 
200 rn 1 1 1 1 1 




-200 1 1 1 1 1 1 

2 4 6 8 10 



FIG. 2 Morse potential V(r) = D[l - exp(-a{r - r e })} 2 - D 
with parameters D = 188.4355, a = 0.711248, r e = 1.9975. 



Eigenvalue plot: Morse case 




-80 - 
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0.1 1 10 100 

Log(IEI) 



FIG. 3 Behavior of the eigenvalue function F(E) with energy 
on a semi-log scale for the Morse potential. 



A. The Multichannel case 

The CFM can be extended to the multichannel (or 
multi-component) problem, where the functions a(E;r) 
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Index 


CFM 


Exact 


1 


-178.798248 


-178.798538 


2 


-160.282181 


-160.283432 


3 


-142.778412 


-142.78006 


4 


-126.287987 


-126.288445 


5 


-110.807388 


-110.808578 


6 


-96.3395233 


-96.3404541 


7 


-82.8832169 


-82.884079 


8 


-70.4389801 


-70.4394531 


9 


-59.0056 


-59.0065727 


10 


-48.5851288 


-48.5854378 


11 


-39.1754532 


-39.1760521 


12 


-30.77771 


-30.7784157 


13 


-23.3919983 


-23.3925247 


14 


-17.0183048 


-17.018383 


15 


-11.6557436 


-11.6559868 


16 


-7.3050122 


-7.30533791 


17 


-3.9661877 


-3.9664371 


18 


-1.6390723 


-1.63928342 


19 


-0.3238727 


-0.32387724 



TABLE II Energy levels of the Morse potential Vir) = 
D[l — exp(— a{r — r E })} 2 — D with parameters D = 
188.4355, a = 0.711248, r e = 1.9975. Middle column values 
are the CFM results whereas the last column values are the 
corresponding exact analytically obtained values. Units are 
such that h = 1 and 2/i = 1. 

and f3(E; r) are no longer scalar functions of E and r. 
Mulichannel case occurs when we have a potential well 
with several local minima or when we are dealing with 
Scattering problems requiring simultaneously incident 
and scattered wavefunctions. In addition, spin leads 
naturally to multicomponent wavefunctions. 

Starting from a Schrodinger system of equations of the 
form: 

JV 

Y,FikVk{r)=Q,i = l i 2,...,N (18) 
fe=i 

where: 

F ik (r) = {2 t i/h 2 )V lk (r) 

F u (r) = {2^/h 2 )(-d 2 /dr 2 -E l + V u (r)) (19) 

The elements Fi k (r) of the matrix F are physical oper- 
ators and the unknown functions y k (r) are partial radial 
waves satisfying the boundary conditions at r — and 

r = oo: 

2/fc(r = 0)=y fc (r = cx))=0Vfc (20) 



Following Friedman and Jamieson (fLsl ) , we rewrite the 
iV-component equation as a set of N coupled Schrodinger 
equations for a matrix potential with elements Vij (r) : 

N 

I/i"(r) + (2^/H 2 )[E - Vu(r)] yi (r) = (2^/H 2 ) £ V lk y k (r) 

fc=i 

(21) 

where i,k = 1, 2, N. 

Using linear superposition, we can write the full solu- 
tion as: 

JV 

y(r)=J2[yk(ro)a ik (E;r)+y' k (roWi k (E;r)} (22) 
k=l 

where, as before, r$ is an arbitrary point between zero 
and infinity, oii k (E; r), f3i k (E; r) are 2N independant par- 
ticular solutions satisfying the set of initial conditions: 

a ik (E;r Q )= 0' ik (E; r ) = 6 ik 

a' lk (E;r ) = f3 ik (E;r ) = (23) 

Going from the indicial to the matrix representation, 
we can build a matrix L(E, r) that depends on the energy 
E and radial position r defined by: 

L(E,r)=/3- 1 (E,r)a(E,r) (24) 

where a and (3 are N x N matrices whose elements are 
the o>i k (E;r),/3i k (E;r) functions. 

The analogue of the eigenvalue function is now given 
by the determinant equation: 

D(E) = |L(E,0) -L(E,oo)| (25) 

The eigenvalues are given by the zeroes of the deter- 
minant D{E) = 0. 

Let us apply the above to the particular case where the 
matrix potential V(r) (fLsh is a constant coupling matrix 
C multiplying a scalar function M(r): 

V(r) = CM(r), Vy = CyM(r), i,j' = l,..6 

The r dependence is contained solely in the radial 
function M(r) = z 2 — 2z, z = exp(— ax), x = r — r e where 
r e is the scalar potential (minimum) equilibrium value. 
C is a full (6x6) constant matrix that may be diago- 
nalised by the similarity transformation: C = QGQ~ X . 
Q is a non-singular real constant matrix and G is a (6x6) 
diagonal matrix. Picking a coupling matrix C similar 
to a diagonal matrix G whose elements Gu,i =1..6 are 
respectively: {24,31.5,37.5,42,56,63}, the initial fully 
coupled six-channel problem [221 is transformed into six 
uncoupled Schrodinger equations that are solved with 
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the multichannel CFM. The full spectrum of V(r) is 
made from the union of the sets of eigenvalues of the six 
individual Schrodinger equations. 

In order to benchmark the CFM results, we note that 
this multichannel example is nothing other than a special 
matrix form of the Morse potential. Then, we have access 
to the full spectrum analytically as in the scalar Morse 
case: 

Ei n = —w e x e (di — n — 1/2) 2 , i = 1, ...6, n < d{ — 1/2, 

(27) 

where di = \J (Gu/w e x e ) and w e x e = a 2 Ti 2 /2/i. In 
order to perform a detailed numerical comparison be- 
tween the numerical performance of the CFM and the 
multichannel Morse potential analytical results, we take: 
r e = X.5A, w e x e = 8 cm -1 and a = 1.540 A. 
The comparison between the computed and the exact 
eigenvalues displayed in table Hill shows that they are al- 
most indistinguishable. 



(i, n) 


Exact 


CFM 


(6,0) 


-42.550 055 679 356 


-42.550 055 679 344 


(5,0) 


-36.833 989 511 483 


-36.833 989 511 484 


(4,0) 


-25.669 697 220 176 


-25.669 697 220 173 


(3,0) 


-22.179 491 924 311 


-22.179 491 924 313 


(2,0) 


-17.625 492 133 612 


-17.625 492 133 614 


(6,1) 


-13.650 167 038 069 


-13.650 167 038 068 


(1.0) 


-12.143 593 539 448 


-12.143 593 539 445 


(5,1) 


-10.501 968 534 449 


-10.501 968 534 446 


(4,1) 


-5.009 091 660 529 


-5.009 091 660 528 


(3,1) 


-3.538 475 772 933 


-3.538 475 772 933 


(2,1) 


-1.876 476 400 837 


-1.876 476 400 834 


(6,2) 


-.750 278 396 781 


-.750 278 396 781 


(1.1) 


-.430 780 618 346 


-.430 780 618 346 


(5,2) 


-.169 947 557 416 


-.169 947 557 416 



TABLE III Analytical and Multi-channel CFM computed 
eigenvalues All values are in cm -1 . 



B. Evaluation of energy spectra and phaseshift 

Before treating any problem with the CFM a number 
of constraints should be underlined in order to properly 
tackle any problem with this method whether it pertains 
to energy spectra evaluation or phaseshift: 

• General considerations: 

— The potential is spherically symmetric. 

— The interval of interest should be of infinite 
length. 



— Ill-conditioning and accuracy: The spectrum 
depends on the zeroes of F(E) — l+(E) — 
l-(E), This subtraction leads to inaccuracies 
because the entire spectrum depends on the 
zeroes of F(E) in the single channel case or 
D(E) = |L(E,0)-L(E,oo)| in the multichan- 
nel case. 

— ro issue and the number of eigenvalues: The 
number of eigenvalues depend strongly on ro- 
ro should be increased until a tan(x) like dia- 
gram for the energy function is obtained. 
For short-range potentials, ro should be close 
to the minimum of the potential or the equi- 
librium value of the radial distance r e . If ro is 
far from minimum, it will see only those levels 
close to the value of V(r e ). 

• Energy spectra calculations: 

— The method being sensitive to convergence for 
r — > and r — > oo, One has to check that for 
r — ► the canonic functions a and f3 are di- 
verging near the origin in the same way (sat- 
uration of the ratio) . 

— One has to check that for r ^ oo the 
canonic functions a and (5 have reached their 
Asymptotic behaviour (of the Coulombic form 
sin(kr)/kr). 

Note: This is not true for the phase. 

— One should have a regular structure of the 
tan(.E) type for the energy function F(E) = 
1+{E) — I -(E) in the single channel (or 
D(E) = |L(E,0)-L(E,oo)| in the multichan- 
nel) case. 

— The method is very sensitive to ro- 

— Forward integration step r — > oo: typically it 
is 0.1 to 1 in energy calculations. In phase 
calculations it is 0.1 to 0.01. It is sensitive to 
the energy used. When the energy increases 
the step should be decreased. There is a com- 
promise to be reached because if the step is 
too much reduced any physical quantity will 
saturate immediately. 

• Phaseshift calculations: 

— The phase stability might be reached long 
before the asymptotic regime in ty(E;r) ~ 
sm(kr)/kr is reached. Typically the phase is 
stable when the potential is about 0.1 au. 

— One has to check that for r — > the canonic 
functions a and (3 are diverging near the origin 
in the same way (saturation of the ratio). 

— The method being very sensitive to ro, we 
ought to treat separately the different types 
of potential. The general rules are enunciated 
below: 



7 



* Short-range potentials: Typically r$ ~ 1 
for all angular momentum L and energies. 
For large energies ro should be reduced. 

* Long-range potentials: Typically ro ~ 1 
for all L but for small energies (1-10 au). 
For large energies r should be reduced to 
0.1. 

— Bessel functions are used for Short-range po- 
tentials only. Trigonometric functions are 
used for Long-range and Short-range poten- 
tials potentials. 

— Forward integration step r — > oo: typically it 
is 0.1 to 1 in energy calculations. In phase 
calculations it is 0.1 to 0.01. It is sensitive to 
the energy used. When the energy increases 
the step should be decreased. There is a com- 
promise to be reached because if the step is 
too much reduced any physical quantity will 
saturate immediately. 

— The stability of the CFM is based on the 
following: two independent sets of solutions 
are generated at some central point, inte- 
grating inwards to the origin and outwards 
to the asymptotic region. Both contain 
linear combinations of the regular and the 
irregular solutions, and by suitably combining 
them, the irregular solution is eliminated. 
Integration can be made out to a very 
large radius, allowing the phaseshift to be 
determined by matching to plane or Coulomb 
wave solutions. It is not necessary to obtain 
a series expansion of the solution in order to 
start the integration, making it unnecessary 
to series expand the potential. 

— The method can be successfully applied to 
scattering problems where the phaseshift is of 
paramount importance in the determination 
of the scattering cross section. Of particular 
importance are (e,2e) scattering problems that 
are dealt with in section VII. 

A final point to consider when dealing with the CFM 
using standard integration methods is that one cannot in 
general proceed toward the origin r = because of the 
potential singularity (for instance, the RK method blows 
up at r = 0) in contrast with the VSCA algorithm. 

III. POTENTIAL ESTIMATION FROM SPECTROSCOPY 

The high accuracy of the CFM allows to retrieve 
the potential (in parametric or pseudo form) from the 
observed atomic or molecular spectra. We exploit the 
Quantum Defect (QD) as an indicator of the results 
obtained by the CFM in order to find the optimal 
parameters of the potential. 



QD information is widely exploited in modern spec- 
troscopy, to characterize Rydberg states and in the 
calculation of the photoionization cross sections of 
various atomic and molecular species (P; [27h . 
We show in this section that QD information, used 
within the Distorted Wave Born Approximation 
(DWBA) framework |62T). might also prove useful for 
the description of ionization processes by presenting 
an alternative way to account for the short range 
interactions (static and exchange) in the calculation of 
the final state continuum distorted waves. The range 
of validity of this approach reaches beyond that ((5lh of 
the commonly used Furness-McCarthy local exchange 
approximation (fl9h . Compared to the determination of 
the Hartree-Fock non-local operator |63h. which becomes 
rapidly a tedious task as the size of the target increases, 
our method allows for a target-independent procedure 
which can be readily applied to much larger atomic or 
molecular systems. 

The parameters of the Green-Sellin-Zachor (|2lh para- 
metric form of the electron-ion potential are optimized 
in order to reproduce the QD using the CFM (|29l). 
Several parametric potentials are discussed extensively 
in the literature and we already used several functional 
forms (0) adapted to different atomic systems. We 
have studied (|59T ) the QD of some rare gases with 
the Klapisch parametric potential ((H) and found 
that in some cases it was very difficult to optimize 
parameters that provide an accurate representation of 
the experimental QD. We believe, the Green-Sellin- 
Zachor (f2lh is better suited to our present study as 
the parameter space is small (two-dimensional) which 
allows for an efficient search of the optimized parameters. 

The electron-ion potentials obtained for each Rydberg 
series are further modified to account classically for the 
electron-electron interaction in the final state. Calcu- 
lations performed within the DWBA framework for the 
ionization of argon in the equal energy sharing geometry 
k a = — kf,, which is reasonably well documented both 
theoretically and experimentally, clearly validate our ap- 
proach and shows significant improvements over previous 
treatments. 



A. Optimisation procedure 

An extensive review of the applications of model po- 
tentials has been given by Hibbert ([24]) and Aymar et al. 
(0). The functional form suggested by Green et al. (f2l[ ) 
is given by: 

V(r) = -(2/r)[{Z-i)u{r) + l], 
with w(r) = l/[e 1 (exp(r/e 2 )-l) + l] (28) 

where t\ and ei are parameters that are determined by 
the optimisation procedure. 
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The optimisation problem, at hand, is over-determined 
since the experimental set of energy levels might consist 
of tens of values whereas the potential depends only on 
two numbers, namely e± and 62 . This over-determination 
allows us to use several criteria for the optimisation 
procedure and later on select the best one that achieves 
results closest to experiment. 

The optimization procedure consists of defining an 
objective function and finding its minimum in the 
two-dimensional parameter space ex, £2- The objective 
function is a quadratic consisting of the difference 
between some picked levels and those produced by the 
parametric potential through the solution of the RSE 
with the CFM. 

We adopted several strategies based on the following 
observations: The QD value is not stable for the low 
levels but tends to reach a stable value when the energy 
increases. Despite the stability of the QD for the 
higher levels, the experimental (and therefore) numerical 
accuracy decreases when the energy increases. 

Therefore a compromise should be achieved by select- 
ing the levels in order to define the objective function to 
minimize. 

We found that a reasonable compromise should be 
based on the following operations that differ with the 
selected rare gas: 

1. Pick some level and take the average of a number 
of higher ones. 

2. Pick two high levels for which the QD has stabilised 
within a given accuracy. 

3. Pick a low level and a high one for which the QD 
has already stabilised. 

All the above operations should yield roughly the 
same values for the parameters before running the final 
check in order to test the accuracy of the obtained 
eigenvalues. 



Series ei 62 
1=0 3.625 1.036 
1=1 3.62 1.06 
1=2 3.6344 1.036 



method for the iterative evaluation of the Jacobian of 
the objective function needed during the minimisation 
procedure. 

It is a Quasi-Newton method that consists of approx- 
imating the Jacobian and updating it with an iterative 
procedure. It converges superlinearly to the solution like 
all secant methods. 

There are several ways to perform the integration of 
the RSE on the basis of the CFM. One may use a fixed 
step scheme such as the explicit RK4 method or a vari- 
able step (VSCA) procedure. Optimization wise, the 
RK4 method is faster than VSCA but less accurate. To 
judge the accuracy of our optimization we compare in 
table [TVl our results for the energy levels and correspond- 
ing QD obtained using the RK4 and VSCA integration 
scheme with those of Szydlik et al. The results dis- 

played in table IIVI clearly favor, as expected, the VSCA 
integration scheme. In contrast to the VSCA, the RK4 
integration scheme is limited to fourth order accuracy. 
Table IIVI further shows the sensitivity of the QD to the 
numerical values of the calculated energy levels. Despite 
small differences between the energy values obtained with 
the RK4 and VSCA methods, the corresponding QD's 
largely differ. This sensitivity of the QD is the main mo- 
tivation for using it rather than the raw energy levels in 
our optimization procedure. 



Experimental levels 


RK4 




VSCA 


-0.309522 


-0.214082 


-0.310563 


-0.124309 


-9.87422(- 


■2) 


-0.124506 


-6.76780(-2) 


-5.69743(- 


■2) 


-6.76904(-2) 


-4.25540(-2) 


-3.70845(- 


■2) 


-4.25546(-2) 


-2.92210(-2) 


-2.60593(- 


■2) 


-2.92238(-2) 


-2.13080(-2) 


-1.93126(- 


-2) 


-2.13062(-2) 


-1.62200(-2) 


- 1.48846 (- 


■2) 


-1.62210(-2) 


-1.27620(-2) 


-1.18223(- 


-2) 


-1.27614(-2) 


-1.03020(-2) 


-9.61654(- 


•3) 


-1.03015(-2) 


-8.49000(-3) 


-7.97521(- 


•3) 


-8.49011(-3) 


-7.11800(-3) 


-6.72092(- 


-3) 


-7.11771(-3) 



TABLE V Comparison between the experimental and cal- 
culated energy levels of the Rydberg series of Argon. The 
levels calculated with the CFM are obtained either with fixed 
step (RK4) or variable step (VSCA) integration. All values 
in Rydbergs. 



TABLE IV Szydlik-Green parameters for the first three Ry- 
dberg series of Argon. 

The optimization program itself is based on a globally 
convergent method for solving non-linear system of 
equations: the multidimensional secant method devel- 
oped by Broyden (|9|). It is based on a fast and accurate 



IV. VIBRATIONAL SPECTRA OF COLD MOLECULES 

A new kind of high precision molecular spectroscopy is 
probing long-range forces between constituent atoms of 
molecules. This spectroscopy is based on using light to 
combine two colliding cold-trapped atoms into a tenuous 
molecule. 
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Experimental QD RK4 VSCA 






202561 





838726 





205576 





163723 





817645 





165967 





156063 





810516 





156415 





152366 





807174 





152400 





150046 





805324 





150326 





149399 





804191 





149110 





148104 





803444 





148345 





148016 





802940 





147808 





147664 





802574 





147425 





147091 





802298 





147161 





147199 





802084 





146957 



TABLE VI Comparison between the experimental and cal- 
culated QD of the Rydberg series of Argon. The levels calcu- 
lated with the CFM are obtained either with fixed step (RK4) 
or variable step (VSCA) integration. 



The burgeoning field of " Photoassociation Spec- 
troscopy" is allowing very precise measurement of life- 
times of the first excited states of Alkaline atoms and 
observation of retardation effects and long-range forces. 
It provides a means of probing accurately the weak in- 
teraction between these atoms (|26l). 

The agreement between theory and experiment re- 
quires simultaneously a highly accurate representation 
of the interaction potential as well as a highly reliable 
method for the calculation of the corresponding energy 
levels. 

Since our aim is directed towards the latter problem, 
we make use of an alternative method to evaluate the en- 
ergy levels for the potential at hand instead of comparing 
to the experimental values in order to assess the validity 
of our results. 

The determination of the vibrational spectra of these 
very tenuous molecules is extremely subtle specially for 
the highest levels which play an important role in pho- 
toassociation spectroscopy. Thus a careful control of ac- 
curacy is needed in order to diagonalise the Hamiltonian 
without losing accuracy for all energies including those 
close to the dissociation limit. 

The magnitudes of potential energy, distance and mass 
values in these kinds of molecules stand several orders of 
magnitude above or below what is encountered in ordi- 
nary short-range molecules. 

For instance, the typical intramolecular potential well 
depth at the equilibrium distance of about 100 clq 
(Bohrs), is a fraction of a cm -1 while the reduced mass 
is several 10,000 electron masses. All these extreme 
values require special numerical techniques in order to 
avoid roundoffs, divergences, numerical instability and 
ill-conditioning during processing. 

As a first example, Johnson et al. l|25h introduced a 
variant of the Morse potential with parameters given by 
D = 0.142, a = 0.815, r e = 2.835. While the mass pa- 



rameter in the Morse case is 0.5 m e , in the Johnson case 
it is very large and equal to 9114.44 m e . 
Large values of the mass parameter lead to a very shal- 
low well in the potential energy (see fig. 0]) that means a 
very small binding energy making the molecule extremely 
weak. This then corresponds to the so-called long-range 
or cold molecule because of the extremely low molecular 
binding energy. In addition the potential width around 
the minimum is very broad implying large quantum fluc- 
tuations in the ground state in contrast to the ordinarily 
deep and narrow Morse potential. 
The potential energy is given by: 

V(r) = D[l - exp(-a{r - r e })} 2 - D (29) 



Johnson potential 



> 




2 4 6 8 10 



FIG. 4 Johnson potential versus distance. The Morse pa- 
rameters are D = 0.142, a = 0.815, r e = 2.835. Note the 
shallowness and breadth of the potential minimum. We use 
the same units as the Morse potential described previously in 

fig.H 

In order to perform the eigenvalue calculation we 
rescale the parameters D = 0.142, a = 0.815, r e = 2.835 
such that they are as close as possible to the Morse poten- 
tial described previously. The scaling yields the param- 
eters: D = 142.387, a = 0.815131, r e = 2.83459. These 
values yield in a straightforward way all the eigenvalues 
(62 in total), the first thirty only being displayed in Ta- 
ble El 

A. Vibrational energy levels of the 23 Na2 molecule 

Since vibrational spectra of the Johnson molecule is 
determined accurately and to arbitrary accuracy with the 
CFM, we move on to tackle the 23 Na2 molecule. The 
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Index 


Exact 


CFM 


Johnson 





497.668182 


497.726562 


498. 


1 


1481.66138 


1481.67773 


1482. 


2 


2449.65503 


2449.69922 


2450. 


3 


3401.64697 


3401.70508 


3402. 


4 


4337.6416 


4337.69531 


4338. 


5 


5257.63672 


5257.68164 


5258. 


6 


6161.63232 


6161.66016 


6162. 


7 


7049.62842 


7049.63672 


7050. 


8 


7921.62256 


7921.67969 


7922. 


9 


8777.62012 


8777.70508 


8778. 


10 


9617.61719 


9617.67578 


9618. 


11 


10441.6133 


10441.6348 


10442. 


12 


11249.6113 


11249.6504 


11250. 


13 


12041.6084 


12041.6914 


12042. 


14 


12817.6074 


12817.6465 


12818. 


15 


13577.6055 


13577.6445 


13578. 


16 


14321.6045 


14321.666 




17 


15049.6045 


15049.6191 




18 


15761.6035 


15761.666 




19 


16457.6035 


16457.6484 




20 


17137.6055 


17137.6523 




21 


17801.6055 


17801.6484 




22 


18449.6074 


18449.6328 




23 


19081.6094 


19081.6562 




24 


19697.6113 


19697.6406 




25 


20297.6152 


20297.666 




26 


20881.6172 


20881.6562 




27 


21449.6211 


21449.6602 




28 


22001.625 


22001.668 




29 


22537.6289 


22537.6328 





TABLE VII First thirty vibrational levels of the Johnson di- 
atomic molecule in cm -1 . Johnson results (limited to the first 
sixteen only) are compared to the CFM and the exact values. 
The corresponding analytical values are almost indistinguish- 
able (see ref. (Hfjl)). The total number of levels is 62 with the 
parameters D = 0.142, a = 0.815, r e = 2.835. 



We start with the g electronic state of the 23 Na2 
molecule. The corresponding potential is given by: 

V(r) = -[(1 - 3X) + v 7 ! - + SIX 2 } (30) 



where X 



9r 3 A ' 



r is the internuclear distance and the parameter C(0 g ) 
is such that: 



lim V(r) 



c 3 (o-) 



(31) 



Identification of the large r limit yields the result: 
^(Og ) — C(0g )/3. In the calculations below, we have 
used 6*3(0") = 6.390 Hartrees.a 3 ) like Stwalley et aL (|57l ). 
The parameter A = 1. 56512. 10 -4 Rydbergs is the atomic 
spin-orbit splitting. Given Cs(0g) and A, the equilib- 
rium internuclear distance is r e — 71.6eto. 

The peculiarities of the 0~ and l u electronic states of 
the 23 Na2 molecule are evident in fig. [5] The l u state is 
higher that the 0~ and the number of vibrational levels 
is smaller because the potential is shallower. This has 
important consequences as we we explore the vibrational 
spectrum below. 




80 120 160 200 240 



energy levels of the g and \ u electronic states of the 
23 Na2 molecule are determined from the Ground state 
up to the continuum limit. The method is validated by 
comparison with previous results obtained by Stwalley 
et al. H3) using the same potential and Trost et al. 
(60) whose work is based on the Lennard- Jones potential 
adapted to long-range molecules. 

The method we use in this work adapts well to this 
extreme situation with the proviso of employing a series 
of isospectral scaling transformations we explain below. 

We apply the CFM to the calculation of the vibrational 
energy levels of a diatomic molecule where the interaction 
between the atoms is given by the Movre and Pichlcr 
potential HI). 



FIG. 5 Potential energy in cm -1 for the 0g and l u electronic 
states (upper shallow curve) of the 23 Na2 molecule. The radial 
distance r is in Bohrs (ao) units. 

We scale all energies with a factor E (usually cm -1 ) 
with the use of equation (1). Then we scale all distances 
with a typical length Lq transforming the RSE appropri- 
ately. This double transformation is reflected generally in 
the potential coefficients preserving thus the functional 
form of the potential. 

In order to gauge the accuracy of the spectra, we per- 
form the integration of the RSE with two different meth- 
ods: the fixed step RK4 and the variable step VSCA 
methods (J35|). 

The VSCA method is based on a series expansion of 
the potential and the corresponding solution to an order 
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-6 1 1 1 1 1 1 1 1 1 1 1 

Ie-10 le-09 le-08 le-07 le-06 le-05 0.0001 0.001 0.01 0.1 1 

Log(IEnergyl) 

FIG. 6 Behavior of the eigenvalue function F(E) with energy 
on a semi-log scale for the l u electronic state of the 23 Na2 
molecule. The vertical lines indicate the eigenvalue position. 
Energies are in cm -1 . 

such that a required tolerance criterion is met. Ideally, 
the series coefficients are determined analytically to any 
order, otherwise loss of accuracy occurs leading quickly 
to numerical uncertainties as discussed later. 

1 I 1 j 1 1 1 1 

0.5 - 




-1.5 - 



-2 - 



-2.5 1 1 1 1 1 1 

40 80 120 160 200 240 

r 



FIG. 7 Potential energy of the Trost et al. HO) model. The 
energy is in e units where e is the potential well depth and the 
radial distance is in where B is a scaled mass. We use the 
same interval for the potential and radial distance in order 
to ease the comparison with fig. \5\ pertaining to the 23 Na2 
molecule. 

Table IVIIII shows the results we obtain with the RK4 
method. The limitation of the RK4 method to fourth 
order hampers the finding of levels beyond the 33rd (see 
table IVlIT|) . In order to find the higher levels we have to 
select an algorithm that enables us to tune the accuracy 
well beyond the fourth order. 

Pushing the accuracy within the framework of a fixed 
step method has the effect of reducing substantially the 



integration step. In order to avoid this problem, we use 
a variable step that adjusts itself to the desired accuracy, 
the VSCA method. 

This method is powerful and flexible enough to find all 
the desired energy levels and allows us to find one addi- 
tional level that was not detected before. It should be 
noted that the last three levels given by Stwalley at al. 
((13) were extrapolated and not calculated. The agree- 
ment between our calculated levels and those of Stwalley 
et al. ((571 ) is quite good. We believe that the small 
discrepancy, increasing as we progress towards the disso- 
ciation limit, is due to a loss of accuracy associated with 
traditional methods in sharp contrast with the CFM. 

The estimation of accuracy of the results hinges basi- 
cally on two operations, integration and determination 
of the zeroes of the eigenvalue function F(E). The 
superiority of the VSCA method is observed in the 
determination of the upper levels that are not detected 
by the RK4 method (see Tables 1 and 2). In addition, it 
is observed in the behavior of the eigenvalue ratio versus 
the index. While in both cases (RK4 and VSCA) the 
ratio increases steadily as the index increases because 
we are probing higher excited states, in the RK4 case 
it rather blows up as dissociation is approached. We 
use typically series expansion to order 12 in VSCA with 
a tolerance of 10~ 8 . In the root search of F(E), the 
tolerance required for a zero to be considered as an 
eigenvalue is 10~ 15 . This does not imply that we disagree 
as strongly as 0.13%, for instance, with the Ground 
state value (see Tables 1 and 2) found by Stwalley et 
al. H3) for the simple reason, we use a splitting energy 
A = 1. 56512. 10~ 4 Rydbergs corresponding to an equi- 
librium internuclear distance r e = 71.6eto. Stwalley et al. 
do not provide explicitly the value of A they use, how- 
ever Jones et al. (|26l ) use a value that is slightly different. 

We move on to treat the 1„ electronic state of the 
23 Na2 molecule. The potential associated with the l u 
electronic state of the 23 Na2 molecule is a lot more in- 
volved. Its VSCA implementation is particularly diffi- 
cult because of the complex functional of the potential 
as we explain below. Analytically, the VSCA algorithm 
requires performing a Taylor series expansion to any or- 
der around an arbitrary point (f3ll ) . 

Numerically, this is still an open problem for arbitrary 
functions. If one turns toward the use of LISP based sym- 
bolic manipulation techniques cumbersome expressions 
are produced hampering any progress. Special methods 
based on analytical fitting expressions are needed in or- 
der to turn the series coefficients into a more manageable 
form (see refs (HI; [13)). 

The first step is to determine the l u electronic state 
of the 23 Na2 molecule by solving the Movre et al. (|42T ) 
secular equation such that: 

V(r) = A[-2 cos(^^) (32) 

C(l ) 

where a = — 2 — 6X and X = n VA ■ In addition, = 
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Index RK4 (cm 1 ) Stwalley et al. (cm 1 ) Ratio 



1 


-1.7864563 


-1.7887 


1.00126 


2 


-1.5595812 


-1.5617 


1.00136 


3 


-1.3546211 


-1.3566 


1.00146 


4 


-1.1704091 


-1.1723 


1.00162 


5 


-1.0057168 


-1.0075 


1.00177 


6 


-0.8592746 


-0.86087 


1.00186 


7 


-0.7297888 


-0.73125 


1.00200 


8 


-0.6159592 


-0.61729 


1.00216 


9 


-0.5164938 


-0.51770 


1.00234 


10 


-0.4301231 


-0.43120 


1.00250 


11 


-0.3556114 


-0.35657 


1.00270 


12 


-0.2917683 


-0.29261 


1.00288 


13 


-0.2374567 


-0.23820 


1.00313 


14 


-0.1916004 


-0.19224 


1.00334 


15 


-0.1531898 


-0.15374 


1.00359 


16 


-0.1212858 


-0.12176 


1.00391 


17 


-9.5022481(-02) 


-9.5438(-02) 


1.00437 


18 


-7.3608067(-02) 


-7.3940(-02) 


1.00451 


19 


-5.6325397(-02) 


-5.6599(-02) 


1.00486 


20 


-4.2530440(-02) 


-4.2754(-02) 


1.00526 


21 


-3.1650256(-02) 


-3.1831(-02) 


1.00571 


22 


-2.3180032(-02) 


-2.3323(-02) 


1.00617 


23 


-1.6679434(-02) 


-1.6791(-02) 


1.00669 


24 


-1.1768423(-02) 


-1.1854(-02) 


1.00727 


25 


-8.1226859(-03) 


-8.1873(-03) 


1.00795 


26 


-5.4687973(-03) 


-5.5165(-03) 


1.00872 


27 


-3.5792655(-03) 


-3.6136(-03) 


1.00959 


28 


-2.2675456(-03) 


-2.2916(-03) 


1.01061 


29 


-1.3831324(-03) 


-1.3995(-03) 


1.01183 


30 


-8.0680818(-04) 


-8.1747(-04) 


1.01321 


31 


-4.461 1287(-04) 


-4.5276(-04) 


1.01490 


32 


-2.3077899(-04) 


-2.3503(-04) 


1.01842 


33 


-9.4777816(-05) 


-1.1252(-04) 


1.18720 


34 




-4.8564(-05) 




35 




-1.8262(-05) 




36 




-5.6648(-06) 




37 




-1.3175(-06) 




38 




-1.9247(-07) 




39 




-1.1215(-08) 




40 




-4.1916(-11) 





TABLE VIII Vibrational levels for the 0J7 electronic state 
of the 23 Na2 molecule as obtained with a fixed step RK4 
method, Stwalley et al. {57j) results and the corresponding 
ratio. Levels 34-40 were not found by the RK4 method due 
to the precision limited to fourth order. 



Index Stwalley et al. (cm x ) VSCA (cm 1 ) Ratio 



1 


-1.7887 


-1.7864488 


1 


00126 


2 


-1.5617 


-1.5595638 


1 


00137 


3 


-1.3566 


-1.3546072 


1 


00147 


4 


-1.1723 


-1.1703990 


1 


00162 


5 


-1.0075 


-1.0057071 


1 


00178 


6 


-0.86087 


-0.8592631 


1 


00187 


7 


-0.73125 


-0.7297908 


1 


00200 


8 


-0.61729 


-0.6159534 


1 


00202 


9 


-0.51770 


-0.5164882 


1 


00235 


10 


-0.43120 


-0.4301217 


1 


00251 


11 


-0.35657 


-0.3556148 


1 


00269 


12 


-0.29261 


-0.2917693 


1 


00288 


13 


-0.23820 


-0.2374560 


1 


00313 


14 


-0.19224 


-0.1916002 


1 


00334 


15 


-0.15374 


-0.1531893 


1 


00359 


16 


-0.12176 


-0.1212854 


1 


00391 


17 


-9.5438(-02) 


-9.5022588(-02) 


1 


00437 


18 


-7.3940(-02) 


-7.3608452(-02) 


1 


00450 


19 


-5.6599(-02) 


-5.6325744(-02) 


1 


00485 


20 


-4.2754(-02) 


-4.2530867(-02) 


1 


00525 


21 


-3.1831(-02) 


-3.1650591(-02) 


1 


00570 


22 


-2.3323(-02) 


-2.3180420(-02) 


1 


00615 


23 


-1.6791(-02) 


-1.6679756(-02) 


1 


00667 


24 


-1.1854(-02) 


-1.1768655(-02) 


1 


00725 


25 


-8.1873(-03) 


-8.1228816(-03) 


1 


00793 


26 


-5.5165(-03) 


-5.4689541 (-03) 


1 


00869 


27 


-3.6136(-03) 


-3.5793742(-03) 


1 


00956 


28 


-2.2916(-03) 


-2.2676168(-03) 


1 


01058 


29 


-1.3995(-03) 


-1.3831806(-03) 


1 


01180 


30 


-8.1747(-04) 


-8.0683859(-04) 


1 


01318 


31 


-4.5276(-04) 


-4.4613249(-04) 


1 


01486 


32 


-2.3503(-04) 


-2.3110168(-04) 


1 


01700 


33 


-1.1252(-04) 


-1.1035443(-04) 


1 


01962 


34 


-4.8564(-05) 


-4.7468345(-05) 


1 


02308 


35 


-1.8262(-05) 


-1.7767388(-05) 


1 


02784 


36 


-5.6648(-06) 


-5.4747950(-06) 


1 


03471 


37 


-1.3175(-06) 


-1.2597092(-06) 


1 


04588 


38 


-1.9247(-07) 


-1.2716754(-07) 


1 


51352 


39 


-1.1215(-08) 








40 


-4.1916(-11) 









TABLE IX Vibrational levels for the 0~ electronic state of 
the 23 Na2 molecule as obtained with Stwalley et al. results, 
the VSCA method and the corresponding ratio. Levels 38, 39 
and 40 of Stwalley et al. |57h are extrapolated with LeRoy 
and Bernstein (|36T ) semi-classical formulae. 



cos x ( 



1+270X-' 
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The parameter C(l u ) is such that: 

C 3 (l t 



lim V(r) 



(33) 
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Index 


RK4 (cm" 1 ) 


Stwalley et al. (cm ] 


Ratio 


1 


0.1319536 


0.13212 


1.00126 


2 


9.0057392(-02) 


9.0192(-02) 


1.00149 


3 


5.9472159(-02) 


5.9574(-02) 


1.00171 


4 


3.7821597(-02) 


3.7896(-02) 


1.00197 


5 


2.3027828(-02) 


2.3080(-02) 


1.00227 


6 


1.3324091(-02) 


1.3359(-02) 


1.00262 


7 


7.2562182(-03) 


7.2787(-03) 


1.00310 


8 


3.6714971 (-03) 


3.6849(-03) 


1.00365 


9 


1.6950002 (-03) 


1.7024(-03) 


1.00437 


10 


6.9533077(-04) 


6.9904(-04) 


1.00533 


11 


2.4102039(-04) 


2.4492(-04) 


1.01618 


12 




6.8430(-05) 




13 




1.3446(-05) 




14 




1.4122(-06) 




15 




3.8739(-08) 




16 




1.2735(-12) 





TABLE X Vibrational levels for the l u electronic state of the 
23 Na2 molecule as obtained with the RK4 method, Stwalley 
et al. results and the corresponding ratio. RK4 found only 
11 levels and levels 15 and 16 of Stwalley et al. are found by 
extrapolation. 



Table IIVI displays the results we obtain with the RK4 
method that cannot find more than 11 levels due to ac- 
curacy limitations. 




lc-07 lc-06 le-05 0.0001 0.001 0.01 0.1 1 
Log(IEnergyl) 

FIG. 8 Behavior of the eigenvalue function F(E) with energy 
on a semi-log scale for the Trost et al. (|60l ) Lennard- Jones 
molecule. The vertical lines indicate the eigenvalue position. 
Energies are in potential depth e units. 



Index Stwalley et al. (cm 1 ~ 


VSCA (cm" 1 ) 


Ratio 


1 


0.13212 


0.13244150 


1.00243 


2 


9.0192(-02) 


9.07598688(-02) 


1.00630 


3 


5.9574(-02) 


6.01645742(-02) 


1.00991 


4 


3.7896(-02) 


3.83982868(-02) 


1.01325 


5 


2.3080(-02) 


2.34584275(-02) 


1.01640 


6 


1.3359(-02) 


1.36187753(-02) 


1.01944 


7 


7.2787(-03) 


7.44243743(-03) 


1.02250 


8 


3.6849(-03) 


3.77999552(-03) 


1.02581 


9 


1.7024(-03) 


1.75281307(-03) 


1.02961 


10 


6.9904(-04) 


7.23013793(-04) 


1.03430 


11 


2.4492(-04) 


2.54856605(-04) 


1.04057 


12 


6.8430(-05) 


7.18249908(-05) 


1.04961 


13 


1.3446(-05) 


1.43120199(-05) 


1.06441 


14 


1.4122(-06) 


1.53931042(-06) 


1.09001 


15 


3.8739(-08) 


4.66022073(-09) 




16 


1.2735(-12) 







TABLE XI Vibrational levels for the l u electronic state of 
the 23 Na2 molecule as obtained by Stwalley et al., the VSCA 
and the corresponding ratio. A new 15th level is obtained 
with the VSCA method. 



Identification of the large r limit yields to the result: 
C 3 (l u ) = C(l u )(v / 7 - 2)/9. We have used in the calcula- 
tions below C3(l u )=1.383 Hartrees.aQ like Stwalley et al. 
The parameter A = 1. 56512. 10 -4 Rydbergs is the same 
as for the 0~ state. 



The next results for the l u electronic state of the 23 Na2 
molecule are obtained with the VSCA method as shown 
in table [V] We find an additional 15-th level in contrast 
to Stwalley et al. who found fourteen and extrapolated 
the last two levels. 

The corresponding graph of the eigenvalue function is 
displayed in Fig. [5] below. 



B. Lennard- Jones molecules 

We apply our methodology to the Lennard-Jones case. 
Our results are compared to the results obtained by Trost 
et al. (60). We start with the levels obtained with the 
RK4 method. The energy unit is e the absolute value 
of potential well depth. The Asymmetric Lennard-Jones 
potential (ALJ) is given by: 



V(r) 



r 



9i 

j, a 



(34) 



It depends on C\ and Ci that yield an equilibrium 
distance at r — r e and a potential depth — e. In order 
to compare with the appropriate litterature, Trost et al. 
(60) use rather the general parameterisation: 



C x = 



[fi -a) 



,0 



[fi -a) 



ftr° 



(35) 



It is scaled in such a way that the energy is expressed 
in units of the potential well depth — e. When a = 6, 



ft 



12 we obtain r e 



and 



4Ci ■ 



and the 



14 



radial distance is in -^2= where B is a reduced scaled 

mass given by B — 2fier 2 . Numerically Trost et al. use 
10 which is the order of magnitude encountered 



B 



in long-range molecules within the framework of their 
system of units. The potential energy in these units is 
displayed in Fig. [7J The RK4 methods yields the results 
displayed in Table EH 



The accuracy limitation of the RK4 results in losing 
the uppermost level 24. Thus we move on to the results 
obtained with the superior VSCA method. All levels are 
obtained with the VSCA and the agreement with Trost 
et al. results is perfect as witnessed by the ratio values 
of Table 6. The eigenvalue function obtained with the 
VSCA method as a function of energy is displayed in 
Fig. 1 



V. TEST OF BOHR CORRESPONDANCE PRINCIPLE 



LeRoy and Bernstein (|36f ) used a WKB approxima- 
tion to derive a semi-classical formula (SCF) for the 
distribution of vibrational levels near the dissociation 
limit of diatomic molecules. If the long-range inter- 
atomic potential is of the form D — C n /r n their formula 
allows the determination of the dissociation energy D 
as well as the tail of the potential (n and C n ) from 
the experimental energies of the highest vibrational 
levels. Furthermore, this analysis opens the possibility 
of testing directly Bohr Correspondence Principle (BCP) 
that states the agreement between semi-classical and 
quantum results improves for larger quantum numbers. 

We consider the rotationless (L=0) 0~ and the l u 
electronic states of the 23 Na2 molecule correlating 
asymptotically to the atomic 2 P^/2 and 2 Sx/2 states. 
The interaction between the atoms is represented by the 
Movre-Pichler |42l ) potential that behaves as — 1/r 3 as 
r — > +oo and we use the same parameters as those of 
Stwalley et al. (f57h . The Asymmetric Lennard- Jones 
potential (ALJ) C\r~^ — Cir~ a is also considered in the 
highly unsymmetric case C2 >> C\ (wc tackle here the 
case (3 — 12 and a = 6, therefore we have asymptotically 
-1/r 6 ) and the Trost et al. © case. Both potentials 
are tailored to the study of Long-range molecules. 

The determination of the vibrational spectra of these 
very tenuous molecules is extremely subtle especially 
for the highest levels which play an important role in 
photoassociation spectroscopy. Special care is needed 
in order to diagonalise the Hamiltonian without losing 
accuracy for all energies including those close to the 
dissociation limit. 

The magnitudes of potential energy, distance and 
mass values in this kind of molecules stand several orders 
of magnitude above or below what is encountered in 
ordinary short-range molecules. For instance the typical 



intramolecular potential well depth at the equilibrium 
distance of about 100 ao (Bohrs), is a fraction of a 
enr 1 while the reduced mass is several 10,000 electron 
masses as observed in the Johnson case previously. 
All these extreme values require special numerical 
techniques in order to avoid roundoffs, divergences, nu- 
merical instability and ill-conditioning during processing. 

Accuracy and its control being of paramount impor- 
tance in this work, the CFM (|29l ) is an excellent can- 
didate because it bypasses the calculation of the eigen- 
f unctions. This avoids losing accuracy associated with 
the numerical calculation specially with rapidly oscillat- 
ing wave functions of highly excited states close to the 
dissociation limit. 

The semi-classical approximation has been already dis- 
cussed in the literature (@; H3). Here we refer to the 
original LeRoy and Bernstein formulation. 

The next section is a discussion of the the Semi- 
classical analysis as formulated by LeRoy and Berstein 
with its implications. We present further the results we 
obtain for the vibrational levels of the 23 Na2 molecule 0~ 
and l u electronic states. We present an additional vali- 
dation of the method with the ALJ molecular potential. 



A. Semi-classical analysis 

The derivation of the LeRoy and Bernstein (36:) SCF 
is based on the WKB condition for the eigenvalues of a 
potential V(r): 



' + 1/2 = 



^E{v) - V{R) dR (36) 



Ri(v) 



where /i is the reduced mass, E{v) is the energy of 
the level indexed by v and Ri(v), R,2(v) are the the 
corresponding classical turning points. 

Near the dissociation limit, one considers E(v) as a 
continuous function of v and by differentiation of eq. 
one obtains: 



dv 



dE(v) 




1 



r i(v) y/E(v) - V(R) 



dR 



(37) 



Using the asymptotic expression of V(r) = D — C n /r n 
and changing the variable of integration to y = R2 (v) / R, 
we obtain: 



dv 



dE(v) 




R 2 {v)/Ri{v) 



y- 2 (y n -iy 1/2 dy (38) 



This integral can be expressed in terms of incomplete 
Euler Beta functions we use to check the accuracy and 
validate the results. LeRoy and Bernstein ((361 ) SCF is 
derived by taking the limit R2{v) / R\{v) — > +00 which 
leads to: 



dE(v) 
dv 



K n [D-E(v)Y 



(39) 
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Index 


RK4 


Trost et al. 


Ratio 


Index 


VSCA 


Trost et al. 


Ratio 


1 


-0.9410450 


-0.9410460 


1.000001 


1 


-0.9410443 


-0.9410460 


1.000002 


2 


-0.8299980 


-0.8300020 


1.000005 


2 


-0.8299963 


-0.8300020 


1.000007 


3 


-0.7276400 


-0.7276457 


1.000008 


3 


-0.7276415 


-0.7276457 


1.000006 


4 


-0.6336860 


-0.6336930 


1.000011 


4 


-0.6336915 


-0.6336930 


1.000002 


5 


-0.5478430 


-0.5478520 


1.000017 


5 


-0.5478480 


-0.5478520 


1.000007 


6 


-0.4698130 


-0.4698229 


1.000021 


6 


-0.4698206 


-0.4698229 


1.000005 


7 


-0.3992870 


-0.3992968 


1.000025 


7 


-0.3992947 


-0.3992968 


1.000005 


8 


-0.3359470 


-0.3359561 


1.000027 


8 


-0.3359533 


-0.3359561 


1.000008 


9 


-0.2794670 


-0.2794734 


1.000023 


9 


-0.2794718 


-0.2794734 


1.000005 


10 


-0.2295070 


-0.2295117 


1.000021 


10 


-0.2295109 


-0.2295117 


1.000003 


11 


-0.1857220 


-0.1857237 


1.000009 


11 


-0.1857222 


-0.1857237 


1.000008 


12 


-0.1477510 


-0.1477514 


1.000003 


12 


-0.1477498 


-0.1477514 


1.000010 


13 


-0.1152270 


-0.1152259 


0.999990 


13 


-0.1152247 


-0.1152259 


1.000010 


14 


-8.776970(-02) 


-8.7766914(-02) 


0.999968 


11 


-8.7766358(-02) 


-8.7766914(-02) 


1.000006 


15 


-6.498640(-02) 


-6.4982730(-02) 


0.999944 


15 


-6.4982534(-02) 


-6.4982730(-02) 


1.000003 


16 


-4.647400(-02) 


-4.6469911 (-02) 


0.999912 


16 


-4.6469838(-02) 


-4.6469911(-02) 


1.000002 


17 


-3.181750(-02) 


-3.1813309(-02) 


0.999868 


17 


-3.1813146(-02) 


-3.1813309(-02) 


1.000005 


18 


-2.059000(-02) 


-2.0586161(-02) 


0.999814 


18 


-2.0585953(-02) 


-2.0586161(-02) 


1.000010 


19 


-1.235370(-02) 


-1.2350373(-02) 


0.999731 


19 


-1.2350173(-02) 


-1.2350373(-02) 


1.000016 


20 


-6.659580(-03) 


-6.6570240(-03) 


0.999616 


20 


-6.6568735(-03) 


-6.6570240(-03) 


1.000023 


21 


-3.048890(-03) 


-3.0471360(-03) 


0.999425 


21 


-3.0470500(-03) 


-3.0471360(-03) 


1.000028 


22 


-1.053690(-03) 


-1.0527480(-03) 


0.999106 


22 


-1.0526883(-03) 


-1.0527480(-03) 


1.000057 


23 


-1.645210(-04) 


-1.9834000 (-04) 


1.205560 


23 


-1.9832 170 (-04) 


-1.9834000(-04) 


1.000092 


24 




-2.6970000(-06) 




24 


-2.6957891(-06) 


-2.6970000(-06) 


1.000449 



TABLE XII Quantum levels of a Lennard- Jones molecule in TABLE XIII Quantum levels of a Lennard- Jones molecule in 

e units, the depth of the potential well as obtained by the e units, the depth of the potential well as obtained by the 

RK4 method, Trost et al. and the corresponding ratio. Only VSCA method, Trost et al. and the corresponding ratio, 
the last 24th level was missed by the RK4 method. 
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FIG. 9 Approach of the semi-classical limit for the QJ and 
l u (upper curve) electronic states of the 23 Na2 molecule with 
the vibrational quantum number. 



where K n is a constant we calculate below for n=3 and 
n=6. Integrating the above equation yields: 

[D - E(v)] 1/2 - 1/n = ^^-K n (v max -v + v D ) (40) 
An 

where v max is the uppermost vibrational level index 
and vd is generally not an integer that depends on inter- 
actions having shorter range (|36l : 1571 ) . 



Taking the origin of energies at the dissociation limit 
and following Gao (|20n. we define a length scale with: 

A. = (^)^ (41) 

This length should be larger than any length scale en- 
countered in the system. In the case n=3, one defines a 
dimensionless bound-state energy by the scaling: 

e s {v) = ~ A 2 E{V) a (42) 
4 ( ft 2 Au )(l//3 3 ) 2 

With the use of eq. (|4"0"f this leads to: 

r / sil/6 [T(l/3)] 

\~ es \ v )\ = 25/331/2^ ( Vmax - v + v d) (43) 

This implies that the spacing between two scaled 
neighbouring vibrational levels taken to power (1/6) is 

a universal constant given by 25/331/2^ numerically equal 
to 1.11292. 

In the Lennard- Jones and the Trost et al. cases (n=6), 
the scaled bound-state energy takes the form: 
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e s (v) 



E(v) 



16 (ft 3 /2/i)(l//%) a 
Proceeding like the previous case, eq. ([4*0|) leads to 

ll/3 _ [r(i/3)] 3 



(44) 



[-es(v)Y 



2 5 / 3 7T 



-(v max - v + v D ) (45) 



In this case, the spacing between two scaled neigh- 
bouring vibrational levels taken to power (1/3) is given 

by 2 5/ 3 7r numerically equal to 1.92763. Thus, one has 
to determine the vibrational levels and examine the limit 
of E(v) as v increases. 



interaction between the atoms is given by the Movre 
and Pichler potential ([26|; [H" 



After scaling all energies with a factor Eq (usually 
cm -1 ) and all distances with a typical length Lq the 
potential coefficients in the RSE are transformed in a 
way such that the functional form of the potential is 
preserved. 



B. Test with the 23 Na 2 molecule 

We apply the CFM to the calculation of the vibra- 
tional energy levels of a diatomic molecule where the 



The eigenvalue function F(E) depicted in fig. [TTil shows 
that the t&n(E) shape is observed again indicating that 
the CFM is able to perform accurately as previously de- 
spite the extreme sensitivity of this problem. 



10 i 1 1 1 1 1 — i — ; rr~i rn — ; — " — ; i; ;;;;;;; — tj 
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Log(|E|) 

FIG. 10 Behavior of the eigenvalue function F(E) with energy on a semi-log scale for the Trost et al. (HO) ALJ (C*i =1,C 2 = 
3000) molecule. The vertical lines indicate the eigenvalue position. Energies are in potential depth e units. 



The integration of the RSE is performed with the potential and the corresponding solution to an order 
VSCA method HH). This method is superior to fixed such that a required tolerance criterion is met. Ideally, 
step methods and is based on a series expansion of the the series coefficients arc determined analytically to any 
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order, otherwise loss of accuracy occurs leading quickly 
to numerical uncertainties. 

This method is powerful and flexible enough to find 
all the sought levels plus additional levels not detected 
before by Stwalley et al. ((571). It should be noted 
that the last three levels given by Stwalley at al. |H3) 
were extrapolated and not calculated. The agreement 
between the rest of our own calculated levels and those 
of Stwalley et al. is quite good. 

We find all sought levels plus an additional level that 
was not detected before. The agreement between our 
calculated levels and those of Stwalley et al. (57) turns 
out to be very good. 

The semi-classical analysis is displayed in Fig. [9] It 
shows clearly that the sequence of Vibrational levels ap- 
proaches the semi-classical limit 1.11292. The departure 
from this limit at the uppermost side of the spectrum is 
stronger for the l u than for the 0~ state. A crossover 
region is spanned by the semi-classical limit and appears 
to be narrower for the l u than for the 0~ case. 

C. Test with Lennard-Jones molecules 

We apply our methodology to the Asymmetric 
Lennard-Jones case (ALJ) as considered by Trost et al. 
(f60h. After obtaining good agreement between the CFM 
eigenvalues and Trost et aZ.'s we display the results of the 
semi-classical analysis in Fig. [6l The sequence of neigh- 
bouring levels approaches the theoretical limit 1.92763. 
The departure from this limit at the uppermost side of 
the spectrum is stronger for the Trost case than for the 
(Ci = 1,C 2 = 3000) ALJ case and the width of the 
crossover region is almost zero in the Trost et al. case. 

Gao (f20h applied this analysis to sequences of neigh- 
bouring levels with his own calculated levels and claimed 
that BCP breaks down for all quantum systems in which 
the asymptotic interaction is of the form 1 jr n with n > 2. 

Additionally Gao found that the agreement between 
his calculated levels and the SCF is better around the 
middle of the vibrational spectra contrary to Stwalley et 
al.'s (H3) and our work. 

Fig. 2 clearly shows that the semi-classical limit is 
reached as we increase the quantum number for the 0~ 
and l u states of the 23 Na 2 molecule. The departure is ob- 
tained with the levels we obtained at the high end of the 
spectra with the VSCA method based on the CFM. The 
crossover region is wider in the 0~ case and the departure 
is smaller. 

On the other hand, Fig. 3 pertaining to the Lennard- 
Jones (V(r) ~ — 1/r 6 ) case shows that the crossover 
region is quite narrow for the ALJ as well as the Trost 
et al. case in comparison with the 23 Na 2 case. The 
departure is also stronger in the Trost case. Additional 



work is needed is order to tie the width of the crossover 
region to the asymptotic character of the potential. 
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FIG. 11 Approach of the semi-classical limit for the ALJ 
(Ci = 1, Ci = 3000) potential as well as the Trost et al. HO) 
model (upper curve) with the vibrational quantum number. 



VI. EXCHANGE: LOCAL AND NON-LOCAL 

Calculations of cross sections for photoionisation re- 
quire the resolution of the integro-differential equations 
satisfied by the radial parts of the free electron wave- 
functions. Very often a simplifying approximation is 
made through replacement of the exchange operator by a 
central equivalent exchange potential (Furness and Mc- 
Carthy ([l9|), Bransden and Noble (0)). While this is 
probably satisfactory for electrons with energies of a few 
eV and more it may cause problems near to threshold. 
Furthermore, to interpret the new generation of exper- 
iments using polarized electrons the proper inclusion of 
exchange may well be important. There exists a means of 
exact solution of the equation for the radial wavefunction 
including static and static- exchange terms (the diagonal 
parts of the direct potential and the exchange operator) 
and a polarisation potential. This is the Distorted Wave 
Polarised Orbital (DWPO) method of McDowell et al 
((39T ) which was developed for Distorted Wave Born cal- 
culations of excitation of Hydrogen (McDowell et al (|40l)) 
and of Helium |56j) and replaces the integro-differential 
equation by coupled differential equations. In both cases, 
only an s-state of the atom is considered and the central 
potentials are expressed in analytical form. 

The phaseshift and the wavefunction in the asymptotic 
region were determined by use of the analytic second or- 
der JWKB solution (Burgess (fiol)'). The solutions are 
started by means of series expansions at the origin and 
continued by Numerov |43l ) integration: the form of se- 
ries corresponding to the regular solution is imposed by 
taking a power r l +1 at the origin in the V partial wave. 
For small k the choice of integration step and changeover 
point is delicate, since convergence of the series requires 
small r whereas the Numerov integration becomes un- 
stable (picking up some of the irregular solution) if it is 
started at too small a radius. 

Accuracy in evaluating the scattering phaseshift 
produced in electron atom collisions is demanded in 
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order to estimate reliably the cross section (j49l ). In (e, 
2e) (|62| ) collision type experiments, triple Differential 
Cross Sections (TDCS) for electron impact ionization of 
an atom or ion, estimated within the Distorted Wave 
Born Approximation (DWBA), are required. 
Computational difficulties such as instabilities, lack of 
accuracy or convergence slow-down arise in solving the 
corresponding Schrodinger equation. The latter has 
to respect Pauli exclusion principle stating that the 
wave function of the incident electron be anti-symmetric 
with the wave functions of the electrons in the target 
atom or molecule. In the Hartree-Fock formulation, this 
requirement leads to the presence of non local terms 
in the resulting coupled differential equations satisfied 
by the radial parts of the free electron wavefunctions 
(0; l62h . The non-local kernel of the integro-differential 
type equation originating from the exchange operator is 
often simplified and a central LEE potential is used or 
some decoupling procedure is employed (|23l). A review 
of the different methods to treat this problem on the 
basis of the Lippman-Schwinger equation is give n by 
Rawitscher et al. (flit). Furness and McCarthy (|l9l) and 
Bransden and Noble (0) cover the DWBA aspects. 

Using a local exchange potential approximation is 
probably satisfactory for electrons with energies of a 
few eV and more, however it may cause problems near 
to threshold. On the other hand, the interpretation of 
experiments using polarized electrons, requires proper 
account for exchange. In addition, it is worthwhile to 
develop a version of the DWBA (|62T ) including non-local 
exchange and applicable to heavy atom targets. This 
stems from the fact, the successful Converged Close 
Coupling (CCC) approach of Bray (@) becomes increas- 
ingly difficult to apply as target complexity increases. 

The Distorted Wave Polarised Orbital (DWPO) 
method is a means for extracting the exact solution of 



the radial wave equation with static, static-exchange 
terms (the diagonal parts of the direct potential and the 
exchange operator) and some polarisati on p otential. It 
has been developed by McDowell et al @ for DWBA 
calculations for the ionization of Hydrogen (McDowell 
et al fiol )) and of Helium fElt). It replaces the integro- 
differential equation by coupled differential equations. 
Only an s-state of the Hydrogen or Helium atoms is 
considered and the central potentials are expressed in 
analytical form. 

The phaseshift and the wavefunction in the asymptotic 
region are determined by use of the analytic second order 
JWKB solution (Burgess (0)). If the energy of the free 
electron is k 2 Ry, the integro-differential equation can be 
rewritten in terms of the variable kr but the value of the 
radius out to which exchange and the short-range part 
of the static potential remain significant, is determined 
by the extent of the electron cloud of the atom. For 
small k, the second order JWKB solution is valid only at 
radii very much larger than that of the asymptotic zone 
where exchange and short-range potentials are negligible. 



A. The integro-differential equation 

We consider the impact of a free electron of energy k 2 
Ry on a one-electron atomic system of nuclear charge 
Z in a Is atomic state of energy Eiq Ry with a ra- 
dial wavefunction Rw(r) = 2Z 3//2 r exp(— Zr). For a 
free electron with angular momentum quantum numbers 
l',m', if we include only on-diagonal potentials and ex- 
change operators and replace the neglected off-diagonal 
coupling potentials by a polarisation potential V po i(r), 
its radial wavefunction Fi/{k,r) satisfies the integro- 
differential equation: 



where: 



and: 



d 2 l'{l' + l) 
dr 2 r 2 



+ k 2 + V ls {r) + V pol {r) + W v (r) 



F,/(fc,r) = 



(46) 



V ls {r) = — --f\R w {r')\ 2 dr' 



*\R w (r')\ 2 dr' 



(47) 



W v {r)F v {k,r) = {-l) s+l R w {r) [E w - k 2 ] 5 V , / R w {r')F v {k,r')dr 



— / R w tf)F v (kS)dT> - / -i?i (r')^(fc/- 

r J0 Jr r 



')dr 



(48) 
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This integro-differential equation can be transformed into the following system of coupled differential equations, 
which is the starting point for McDowell et al (|39t |4(J): 



+ (-l) S i? 10 (r) 



dr 2 
2 1 



-Vu(r) 



F v (k,r) 



r 21' + 1 



G v {k,r) = {-l) b+1 R w {r)8 Vfi A{k) 



3 2 1/(1' + 1) 21' + 1 

jj^G v (k, r) - K £ } G V (k, r) + -^—R w {r)F v (k, r) = 



(49) 



r 



where 

/>oo 

A(k) = [k 2 - E w ] / R 10 (r')F l ,(k,r')dr' (50) 
Jo 

is a term appearing for the special case of zero orbital 
momentum (/' = 0). Si^q is the Kronecker delta and 



Vi,(r) = k 2 + V ls (r) + V pol (r) 



(51) 



Our aim is to solve this system subject to the boundary 
conditions |39t |40| ) for the Hartree functions: 



Fj,(fc,r) ; 
F r (fc,r) 



Gvik.r) 



a v {k){si.{kr)-ta,n[5ir{k)]ci.{kr)} G v (k,r) — » 

j r — >-oo 



(52) 



where a/'(fc) is a normalisation factor, <5;'(fc) is the 
phaseshift for specific {fc,^} and s;/(/o) and c;<(p) are re- 
spectively: 

• p multiplied spherical Bessel and Neumann func- 



• regular and irregular Coulomb wavefunctions when 
Z > 1 (so that Vi s (r) is a long range potential 
behaving like (Z — l)r _1 when r tends to infinity). 



In the more general case of an ion with a frozen core 
tions when Z = 1 (so that Vi s (r) is a short range and an outer shell of electrons in {n, 1} states of ra- 
potential falling off faster than as r tends to dial wavefunction R n i(r), and a free electron in the state 
infinity); {fc, I'} we have a larger set of coupled equations: 



where: 



dr 2 



F v {k,r) 



V n i{r) 



F v {k,r) 



-(-l) s R nl (r)-y2ju',xGt(k,r) = {-l) S+1 R nl (r)5 LV A nir {k) 
x 

|1 G a r) (fc, r) + ^-Rm (T)F V (M = 



(53) 



V n i{r) = k 2 -V nl (r) 



{r r r 'A poo A ~i 

7 o ^TT\ R ni(r')\ 2 dr' + J \R nl (r')\ 2 dr' j + V pol (r) 



(54) 
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and: 



AniM(k) = [k 2 - E nl ] / R^Fvik^dr' 



These equations are subject to the boundary conditions, for all possible A values: 



F v (k,r)—>Q 

r— >0 



(55) 



GHKr)— >0 

r—*0 



F v (fc, r) — > (k) {s v (At) - tan ft, (fc)] cj/ (fcr)} (fc, r) — > 

r — udo r— »oo 



(56) 



^/.c,A and Ji i> \ are angular integrals which depend on 
the number of electrons in the ion outer shell and the 
angular momentum coupling scheme. V n i (r) is a central 
potential for attraction of an electron by the core and E n i 
is the total energy of the outer shell electrons. This is not 
applicable to hydrogenic ions as the degeneracy of the 
energy in / makes it essential to include channel coupling 
potentials. We will not consider it further here except to 
note that Ji t i' t \ imposes the triangular rule \l — l'\ < A < 
I + I' and I + 1' + A even. This gives an idea of the number 
of different Gf, present and the extent of the problem 
we ultimately wish to solve, in the case of more complex 
atoms. 



B. The CFM for solving the DWPO equations 

In order to facilitate the presentation it is convenient 
to use /i in place of F ; < and f 2 in place of GV so as to 
rewrite the coupled equation system (2) as a special case 
of the more general system: 



f{{r) + V 11 (r)f 1 (r) +V 12 {r)f 2 {r) 
f2(r) + V 22 (r)f 2 (r)+V 21 (r)f 1 (r) 



= 5 l ,, A(k)W 1 (r) 
= S v , A{k)Wa(r) 

(57) 



with 



J 



Vu{r) = V u {r) - 1 ^ V 12 (r) = {-l) s R w {r) 



V 21 (r) 



(21' + 1) 



i?ioW, V 22 (r) 



r (21' + 1) 



Wi(r) = (~l) b+i R w (r), W 2 (r) = 

We can construct the general solution of the above equations as (see Kobeissi and Fakhreddine H2)): 

fi(r) = /i(r )an(r) + f[(r )(3 n (r) + f 2 (r )a 12 (r) + ^(r )/3 12 (r) + S t > fi A(k)at(r) 
h(r) = h(r )a 2l (r) + f{(r )(3 21 (r) + f 2 (r )a 22 (r) + f^(r )(3 22 (r) + S v , A(k)a 2 (r) 



(58) 



(59) 



r 



where {a\j(r), ot 2 j(r)} and {f3\j(r), (3 2 j (r)} are two dif- 
ferent pairs of independent solutions of the homogeneous 
system: 

g'ir) + Vi 1 (r)gi(r)+V a {r)g 2 {r) =0 (60) 
g' 2 \r) + V 22 (r)g 2 (r) + V 21 (r) gi (r) =0 (61) 

satisfying the initial conditions at an arbitrary point r = 

oaj (ro) = Plj (ro) = St j , aty (r ) = fcj (r ) = 

and {tJi(r), cr 2 (r)} is a particular solution of the inhomo- 
geneous system: 



h' 2 '(r) +V 2 2(r)hz(r) +V 2 x(r)hx(r) = W 2 (r) (63) 
satisfying the initial conditions at r = ro 



a i (r )=<7' i (r ) = (64) 



The general solution in matrix form is written as: 



Y(r) = a(r)Y(r ) + P(r)Y'(r ) + 5 Vfl A(k)a(r) (65) 



K(r) +V u (r)fti(r) + V 12 (r)h 2 (r) = Wi(r) (62) with 
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Y(r) 




, a{r) 



otu(r) a 12 (r) 
»2i(r) a 22 (r) 



/3n(r) /?i 2 (r) 
/?2i(r) /3 22 (r) 



(66) 



r 



Each column of the (2x2) matrices a(r) and (3(r) 
is a particular solution of [HI] with initial values given 
by 1621 The column matrix:er(r) = ( a\ (r) , a<z (r) ) is 
a particular solution of [53] with initial condition given 
by [64] Using the boundary conditions [56] imposes: 
a(0)y(r ) + 0(O)Y'(r o ) + 8 V A(k)a(0) = 0, which leads 
to /3- 1 (0)a(0)y(r ) + F'(r j + ^W/^OMO) = 
where /3 _1 (r) is the inverse of the matrix f3(r). 

Thus, the constant matrices at point r = r$, Y'(ro) 
and Y(ro) are related by: 



Y'(r )=Y(r )A + S l , fi A{k)\, 
where the matrices A and A are given by: 



(67) 



The initial values at the arbitrarily chosen starting 
point r = tq are: 

cp(r ) = I, <p'(ro) = A, 7 (r ) = 0, 7 '(r ) = A (71) 



where / is the unit matrix. The solution constructed 
from ip(r) and j(r) is a particular solution of the coupled 
equations (8) for which the functions {/i(r), / 2 (r)} are 
regular at the origin. 



A = -/3- 1 (0)a(0), A = (0)ct(0) (68) 

Substituting back into [65] we then get: 

Y(r) = tp(r)Y(r Q ) + S Vfi A(k)j(r) (69) 

the functions ty?(r) and 7 (r) are related to the particu- 
lar solutions a(r),(3(r) and a(r) by: 

<p(r) = a(r) + /3(r)A, 7 (r) = /3(r)X + a(r) (70) 

We notice again that ip(r) — ( 0i(r), 2 (r) ) and 7 (r) — 
( 7i ( r ) j 72 ( r ) ) ar e particular solutions of system (8) since 
they are linear combinations of a(r) and (3{r). 



C. Phaseshift calculation 

From the first of boundary conditions we can determine 
the phaseshift Si from: 



Q(r) 



tan 5;/ = lim Q(r), where: 

r— *oo 

mr)s,.(kr)-f 1 (r)ks' l ,(kr) 
~ f{{r)c v {kr) - fx{r)kc' v {kr) 



(72) 



We follow Kobeissi et al (|35l ) in using the recursion 
relations for sy (p) and cy (p) : 



Q(r) 



Z = 1 





fi(r) 


si'(kr) - 




S;/ + i(fcr) 


[m ^ 


-fi(r) 


Q'(fcr) - 


hfc/i(r 


e/' + i(Ax) 



Z > 1 



Qir) 







s v (kr) + 


\A 2 + f^ 




su + i(kr) 






cp{kr) + 


^ 2 + {^ 




Q' + i(fcr) 



(73) 



The function Q(r) can be obtained for any radius r. 
We calculate it at large r values and examine its be- 
haviour. When it tends to a constant limit, we con- 
sider that the asymptotic region has been reached and 
the phaseshift is determined modulo 2ir. 

Numerically, the phaseshift is actually calculated with 



the following steps: 

1. An arbitrary point tq is chosen as the origin to 
proceed with the integration. 

2. The canonical functions a(r), f3(r) and tr(r) are 
computed for r < r$ : Eqs l61l and l63l are integrated 
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starting at with a(r), (3{r) and a(r) initialised 
with [^21 and [Ml The matrices A and A are deduced. 
Progressing toward the origin r — 0, the integration 
is stopped when the matrices A and A attain stable 
values. Numerically, this is done by monitoring the 
values of their determinants. 

3. The computation of A and A serves to calculate 
the other canonical functions <p(r) and j(r) and 
consequently estimate A\ and Ai. 

4. Repeating the computation with r > tq allows to 
calculate the ancillary functions fi(r), f[(r) and 
Q(r) (relations [72] and 175)) . The calculation is 
stopped when Q(r) attains a stable value corre- 
sponding to tan 5. The stability of the phase- 
shift value is attained when the ratio |Q(?V|-i) — 
Q(r p )\/\Q(r p+1 ) +Q(r p )\ < e s , where [r p ,r p+ i] is 
the interval obtained after p integration steps and 
eg is termed the stability error. 

The phaseshift overall error (desired error or tolerance) 
depends on both ex and es & s well as computer/compiler 
arithmetic, roundoff errors and intermediate algorith- 
mic/numerical operations. We estimate it through a 
comparison with analytical or other available results 
(such as asymptotic or at special points). 

For small k the outward Numerov integration of the 
regular solution can get out of control if it picks up 
a tiny fraction of the irregular solution, because of ill- 



conditioning due to round-off errors. An alternative is to 
modify the potential by introducing a hard core: V„z(r) 
is set artificially to infinity for r < r s , where r s is the 
starting point of the integration and retains its original 
form for r > r s . 

As mentioned by Bayliss (4) et al, this method gives re- 
sults sensitive to the point in the classically forbidden 
region at which the integration is started. If the starting 
point is too small some solutions become unstable; if it is 
too large for the initial conditions used, the solution is in- 
accurate. In the CFM, the solutions a(r), (3(r) and a(r) 
initially generated are in fact linear combinations of the 
regular and irregular solutions of the coupled equations 
and by taking the linear combinations tp(r) and 7(7*) we 
eliminate from them the irregular solution. 



D. Phaseshift accuracy and comparison to the S-IEM 

We apply the present method to the case of the colli- 
sion of a low-energy electron with atomic hydrogen, i.e. 
generating the wavefunctions needed for DWBA calcula- 
tions of electron impact ionization of atomic hydrogen. 
In this case the static potential is given by: 



V ls (r) =-2 1 



1 



exp(— 2r) 



(74) 



Like McDowell and collaborators, we use a Callaway- 
Temkin polarisation potential (see Drachman and 
Temkin ([12))) of the form: 



J 



9 



1 - e- 2r I 1 



2r + 2r z 



—r 
3 



— r 
3 
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(75) 



r 



To integrate the coupled equation system, preference is 
given in the present work to the VSCA method described 
above. The latter method was shown by Kobeissi et al 
(fill l3Cl l33l ) to be highly accurate in the case of both 
single and coupled differential equations. It requires po- 
tentials to be expressed in analytical form. Numerical 
potentials can generally be fitted by piecewise analytical 
or polynomial functions such as cubic splines. The in- 
tegration can be safely taken out to a very large radius, 



where Vi s (r) assumes its asymptotic form and the phase- 
shift determined. We refer to this method as the CFM 
Exact Exchange (CFMEE). 

In order to further gauge the reliability of our results, 
we performed a series of calculations varying the stability 
error es and the truncation error ex in order to test the 
accuracy and convergence speed of our method. The tests 
are displayed in the tables below for the singlet case: 



r 



Recently, Rawitscher et al. (j49l ) obtained very accu- 
rate results for the phaseshift at I' = for momentum 
k = 0.2/ ao in the absence of any polarization potential 
but with rigorous inclusion of the Fock exchange term. 
They obtained the phaseshifts and scattering lengths in 
the singlet and triplet states. The method they use, the 



S-IEM (Spectral Integral Equation Method), is based on 
partitioning the integration interval [0, r max ] into a fixed 
number of partitions expressing the integration kernel 
with a set of four different functions. This transforms 
the problem into a block tridiagonal matrix system built 
from sets of 4x4 block matrices corresponding each to a 
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TABLE XIV Singlet (S = 0,1' = 0) phaseshifts for e-H scattering versus momentum as the stability error is decreased while 
truncation error er is fixed at 10~ 14 (upper sub-table) and as the truncation error is decreased with a fixed stability error 
eg = 10 -14 (lower sub-table) 



fc 

0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
1.0 



fc 

0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
1.0 



Stability error es 



10" 



10 



-6 



lo- 



ur 



10" 



10" 



2.54421122 
2.06575195 
1.70757716 
1.42957123 
1.2095031 
1.03735556 
0.90616322 
0.80734542 
0.75274001 
0.69239713 



2.53027301 
2.03935101 
1.67206099 

1.3914665 
1.17365026 
1.00552357 
0.87811171 
0.78365856 
0.75420532 

0.7010382 



2.5274412 
2.03407088 
1.66518977 
1.38497592 
1.16825708 
1.00072357 
0.87375939 
0.77961723 
0.71342049 

0.6701203 



2.52744125 
2.03407088 
1.66518968 
1.38497592 
1.16825708 
1.00072404 
0.87376024 
0.77961723 
0.71342026 
0.6701203 



2.52744125 
2.03407088 
1.66518968 
1.38497592 
1.16825708 
1.00072404 
0.87376024 
0.77961723 
0.71342027 
0.6701203 



2.52744125 
2.03407088 
1.66518968 
1.38497592 
1.16825708 
1.00072404 
0.87376024 
0.77961723 
0.71342027 
0.6701203 



Truncation error er 



10" 



10" 



10" 



10 



-10 



10 



-12 



10" 



3.00127274 
2.82900282 
2.65292424 

2.4734271 
1.15828993 
2.14501974 
0.86558783 

1.8745865 
1.72900321 
1.60803611 



3.0307537 
2.88186498 
1.65536145 
1.37484701 
1.15870411 
0.9920774 
0.8658675 
0.77331792 
0.70870253 
0.66656396 



2.52412523 
2.02681296 
1.65611743 
1.37567631 
1.15948866 
0.99279684 
0.86653892 
0.7738408 
0.70910062 
0.66686498 



2.52463171 
2.02795407 
1.65755527 
1.37717447 
1.16092162 
0.9941118 
0.86774235 
0.77481281 
0.70982718 
0.66741118 



2.52558189 
2.03007799 
1.66025047 
1.37993599 
1.16354419 
0.99649203 
0.86992432 
0.77655005 
0.71113237 
0.6684008 



2.5274412 
2.03407088 
1.66518977 
1.38497592 
1.16825708 
1.00072357 
0.87375939 
0.77961723 
0.71342049 

0.6701203 



given partition |49h. 

The phaseshift expressed in the form fc/ tan(<5) versus mo- 
mentum fc is displayed for the singlet and the triplet in 
Fig. 5. Our results (crosses) fall exactly on top of the 
continuous curves for the singlet and the triplet obtained 
by Rawitschcr et al. © with the S-IEM. Whenever, ei- 
ther of the curves cross the horizontal axis, the phaseshift 



value is 7r/2 modulo (it) and tan(<5) is singular. The CFM 
method like the S-IEM, had no particular difficulty (in- 
stability or slowing down) in reproducing the singularity 
and the phaseshifts were obtained as described in sec- 
tion 5.1, for arbitrarily small values of the momentum k 
(k > 1CT 8 in l/oo units). 



The rate of convergence as a function of stability error 
es of the CFM result for the triplet phaseshift at fc = 
0.2/ao is displayed in table IXVl along with the S-IEM 
result. It is remarkable to observe that as the stability 
error is decreased, we obtain steady convergence toward 
the S-IEM result. If we consider, in that case, the S- 
IEM result as exact, we can infer that the overall desired 
error (tolerance) behaves as ~ es ' 45 , implying that es 
is a good indicator of accuracy despite its deterioration 
across intermediate arithmetic operations. 



The scattering length a for the electron-hydrogen atom 
scattering is based on the low momentum expansion of 
the scattering phaseshift: 



k cot S = -l/a + r e k 2 + 0(k 3 ) (76) 

where r e the effective range (j49| ) . The left-hand side of 
the above expression is calculated for two very small val- 
ues of the momentum fc differing by a factor of two and 
the values of a and r e are solved for pH). This procedure 
is repeated for decreasing values of fc and stability error 
£5 until the numerical value of the result is obtained to 
a given number of significant figures. In table IXVI1 we 
display the results of the CFM and S-IEM for the sin- 
glet case. The agreement between the S-IEM result for 
the scattering length 8.100312397 and the corresponding 
CFM result (8.100307932) is excellent. 
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FIG. 12 Singlet and Triplet I' — 0, CFMEE results (crosses) values of fc/tan(<5) versus incident momentum k in 1/ao units and 
comparison to Rawitscher et al. (continuous lines). The smallest value used for k is 10 -8 . 



TABLE XV Triplet phaseshift accuracy as the stability er- 
ror is decreased at k — 0.2/ao. S-IEM is Rawitscher et al. 
Spectral Integral Equation Method 

Stability error es S 



10" 


-4 


1.8927821 


10" 


-5 


1.8773093 


10" 


-6 


1.8724465 


10" 


-7 


1.8708870 


to- 


-8 


1.8703880 


ur 


-9 


1.8702301 


10 _ 


12 


1.8701771 


S-IEM 


1.8701579 



TABLE XVI Singlet scattering length versus momentum as 
the stability error is decreased. S-IEM is Rawitscher et al. 
Spectral Integral Equation Method 







Stability error 


es 


k 




io- 4 


KT 6 


IO -8 


10" 


2 


8.058030140 


3. 104593677 


8.109676748 


lO - 


3 


8.081922310 


3.098359718 


8.100205405 


io- 


4 


8.093862531 


3.099570349 


8.100235983 


10" 


5 


8.098414768 


3.100045256 


8.100281954 


io- 


6 


8.100073345 


3.100226775 


8.100299303 


io- 


8 


8.100305162 


3.100305163 


8.100307932 


S-IEM 






8.100312397 



The CFM result is obtained for a value of momentum 
k = 10~ 8 /ao and stability error es — 1CP 8 , whereas the 
S-IEM result was obtained with a value of momentum 
k = 10~ 5 /ao, and maximum integration range r max ~ 50 
with approximately 1000 mesh points. 
Finally we display in table IXVLTl the result for the triplet 
scattering length as momentum k decreases steadily to- 
ward k — 10~ 8 /ao. 



TABLE XVII Triplet scattering length (with 10 -8 stability 
error) as momentum decreases and comparison to the S-IEM 
result. Note the faster convergence than the singlet case. 



k 




a 


1(T 


2 


2.3497331 


io- 


3 


2.3493995 


io- 


4 


2.3493961 


10" 


5 


2.3493961 


io- 


(i 


2.3493961 


io- 


7 


2.3493961 


10- 


8 


2.3493961 


S-IEM 2.3493961 



Roche et al. I|5a l developed, recently, an exact ex- 
change method, within the DWBA, for the evaluation 
of high and low energy electron impact ionisation ex- 
periments on Hydrogen and Lithium, with presence of 
polarisation. In the table below, we compare our re- 
sults (CFMEE) to theirs for the singlet and triplet s-wave 
phaseshifts at high energies. The polarisation used is of 
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the form: V po i(r) 



VII. ONE-DIMENSIONAL POTENTIAL PROBLEMS 



TABLE XVIII Singlet s-wave phaseshifts with polarization 
V po i(r) — — o 2 ! and comparison to the CFMEE method. 



Energy in eV 


CFMEE 


Roche et al. (53) 


14.6 


-1.46166 


-1.46084 


15.0 


-1.4720848 


-1.47127 


15.6 


-1.48717379 


-1.48642 


17.6 


-1.53373849 


-1.53297 


20.0 


- 1.58264557 


-1.51896 


25.0 


-1.66616627 


-1.66483 


30.0 


-1.73100847 


-1.73061 



TABLE XIX Triplet s-wave phaseshifts with polarization 
V po i(r) = — an d comparison to the CFMEE method. 

Energy in eV CFMEE Roche et al. (53) 



14.6 


-1.45993029 


- 1.43966 


15.0 


-1.46697595 


- 1.44676 


15.6 


-1.47691298 


- 1.45683 


17.6 


-1.50572363 


- 1.48681 


20.0 


-1.53328328 


- 1.51529 


25.0 


-1.57639125 


- 1.56060 


30.0 


-1.60975092 


- 1.59548 



The remarkable agreement shows that our method per- 
forms well at high energies and presence of polarisation, 
with methods using exact exchange instead of local ex- 
change. Our work has also shown that a distorted wave 
code, with exact treatment of exchange and Numerov 
integration, breaks down with appearance of numerical 
instabilities at energies below about 5 eV. Codes with 
LEE potentiels give poor results for low partial waves at 
energies below 15 eV. 
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In the last section before the conclusion, we apply the 
CFM to a couple of ID potentials to illustrate once again 
its capabilities beyond the focus of this work mainly con- 
cerned with the 3D case. 

The archetypical ID problem is the infinite square well 
potential of width w defined by: 

V(x) = 0, < x < w, V(0) = oo, V(w) = oo, 



the eigenvalues are given by: E„ 



tl t 717T \ 

2m„ V w ' 



1,2. 



with m p is the particle mass (not necessarily an electron 
as generally assumed). 

In order to conform with the CFM general framework, 
one transforms the closed finite interval [0, w] to the open 
infinite interval [0, oof. 

It is interesting to note that while, in general, the CFM 
maps an arbitrary potential (infinite interval) over the 
infinite square well potential (finite interval), this prob- 
lem might be viewed exactly as the reverse mapping. 
With the use of the transformation r = ,„ x ^ the 
Schrodinger equation is transformed into: 



(w- 



h 2 d 2 ^(x) 
2m„ dx 2 



Eip(x), < x < w 



(77) 



with the boundary conditions: ^(0) = 0,ip(w) = 0. 
The Schrodinger equation defined over the interval < 
r < oo acquires a first-order derivative term becoming: 



2m Ti 



(1 + r) A d 2 ^(r) 2(1 + rf <#(r) 



dr 2 



dr 



Ei>(r), 
(78) 

The VSCA algorithm derived earlier for the standard 
Schrodinger equation must be altered into (using units 
such that h = 1,to p = 1): 



(n + 2)(n + l)C% = -2 £ (m + l)a£2 m C&>_ i - 2Ew 2 £ 6 £> ro C$ri 

I 



(79) 



m=0 



m— 



where the cim\bm^ are, respectively, the coefficients of 
the series expansion of ■ 



and 



over the interval 



(l+r) (1+r) 4 

r P ,r p+ i]. 

The results are displayed in the table below proving once 
again the remarkable versatility of the CFM method. 

After proving the success of the CFM over a finite 
interval we move on to the standard application over the 
usual infinite interval [0, oof. 



Double-minimum Potential Well (DPW) problems 
are interesting to solve as they arise in many areas 
of Atomic, Molecular and even Solid State physics. 
They are encountered in highly excited Rydberg states 
of atoms in crossed electric and magnetic fields and 
in certain molecular potential curves. When two- 
dimensional electron layers (such as in semiconducting 
hctcrostructures) are placed in perpendicular electric 
and magnetic fields, a potential well with two minima, 
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CFM Index Exact Ratio 






.004 


440 


783 


1 





.004 


444 


444 





.9992 





.017 


917 


823 


2 





.017 


777 


776 


1 


.0078 





.040 


324 


728 


3 





.039 


999 


999 


1 


.0081 





.071 


671 


344 


4 





.071 


111 


105 


1 


.0078 





111 


964 


293 


5 





.111 


111 


104 


1 


.0076 





.161 


208 


108 


6 





.159 


999 


996 


1 


.0075 





.219 


427 


973 


7 





.217 


777 


774 


1 


.0075 





.286 


566 


615 


8 





.284 


444 


422 


1 


.0074 





362 


612 


218 


9 





.359 


999 


985 


1 


.0072 





.447 


584 


212 


10 





.444 


444 


418 


1 


.0070 





.541 


484 


714 


11 





.537 


777 


781 


1 


.0068 





.644 


318 


163 


12 





.639 


999 


986 


1 


.0067 





.756 


076 


396 


13 





.751 


111 


090 


1 


.0066 





.876 


762 


569 


14 





.871 


111 


095 


1 


.0064 


1.006 386 52 


15 


0.999 


999 


| 94 


1 


.0063 



TABLE XX Quantum levels of the infinite square well poten- 
tial given by the CFM along with exact results and the ratio 
of the CFM to the exact eigenvalue. The well width w is cho- 
sen in a way such that the eigenvalue is 1 when the number 
of levels is 15. 



for the electronic motion normal to the surface, arises. 

A DPW can be symmetric or asymmetric and one has 
to adapt in each case the appropriate boundary condition 
imposed by the CFM (see Appendix on matching and 
boundary conditions). 

An interesting DPW is the symmetric Double Gaussian 
potential investigated by Hamilton and Light (22). It is 
given by: 

V(r) = -£)[exp(-fi(r - r a f) + exp(~Q(r + r a ) 2 )] 

The values of the parameters: D, 51, r a are respectively: 
12.0,0.1,5.0 in standard a.u such that K = 1, n = m e = 1. 

The results shown for all the 24 levels in table IXXII 
prove unambigously the accuracy of the CFM method 
and its high level of competitivity with respect to the 
somehow sophisticated method used by Hamilton and 
Light ((22T ) based on Distributed Gaussian Basis sets in- 
spired from Quantum Chemistry Techniques. 

Johnson has studied an asymmetric DWP consisting 
of the sum of a Morse and a Gaussian resulting in the 
expression: 

V(r) = D[l - exp(-B(r - r a ))] 2 + Acxp {-C(r - r h f) 

The values of the parameters A, B, C, D, r a ,ri, are (fol- 
lowing Johnson |25T)) in (cm -1 , A system of units) are: 
10 4 cm -1 , 1.54 A" 1 , 200.0 A" 2 , 31250.0 cm" 1 , 1.5 A, 1.6 
A respectively. 

Eigenvalues for the asymmetric double minimum po- 
tential problem are given in table IXXIII and a compar- 



CFM Index Hamilton and Light 



-11.250 418 469 





-11.245 199 313 


-9.779 202 


594 


2 


-9.773 


496 


902 


-8.387 701 


732 


4 


-8.381 


307 


510 


-7.079 415 


041 


6 


-7.072 


039 


562 


-5.858 805 


474 


8 


-5.849 


958 


02 


-4.732 231 


858 


10 


-4.720 


829 


36 


-3.709 907 


559 


12 


-3.694 


518 


38 


-2.807 436 


691 


14 


-2.798 


251 


92 


-2.022 064 


904 


16 


-2.089 


661 


3 


-1.293 090 


067 


18 


-1.462 


202 


9 


-0.601 483 


056 


20 


-0.771 


081 




-0.067 153 


689 


22 


-0.177 


181 




-11.250 421 409 


1 


-11.245 199 313 


-9.779 225 


834 


3 


-9.773 


496 


902 


-8.387 719 


137 


5 


-8.381 


307 


491 


-7.079 412 


929 


7 


-7.072 


038 


846 


-5.858 811 


221 


9 


-5.849 


940 





-4.732 171 


001 


11 


-4.720 


509 


6 


-3.709 113 


861 


13 


-3.690 


475 


6 


-2.801 628 


760 


15 


-2.763 


219 


7 


-2.000 566 


637 


17 


-1.924 


577 




-1.255 332 


005 


19 


-1.149 


254 




-0.561 216 


170 


21 


-0.457 


88 




-0.045 810 


537 


23 


-0.003 


41 





TABLE XXI Computed eigenvalues for the symmetric dou- 
ble Gaussian well potential. The even index results are shown 
separately from the odd index since they correspond to differ- 
ent symmetries. The numbers at left are the levels computed 
with the CFM^on the right the levels obtained by Hamil- 
ton and Light l|22i ). Note the deterioration of accuracy of the 
Hamilton and Light results as the index increases because of 
the approach of the continuum. 

ison between Johnson's ((25m results and the CFM are 
displayed below. 



VIII. CONCLUSIONS 

In this review, we have shown that the CFM is a very 
powerful and accurate method that is able to solve a large 
variety of quantum problems such as: 

• Energy levels for regular and singular potentials 
(radial case) 

• Pseudo-potential estimation (parametric potential) 
from spectroscopic data 

• Accurate Phase shift evaluation for regular and sin- 
gular potentials 

• Test of the Bohr Correspondence Principle 
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Index 


Johnson 


CFM 





1302.500 


1302.498 972 


1 


3205.307 


3205.303 782 


2 


4227.339 


4227.336 543 


3 


5144.251 


5144.243 754 


4 


6064.241 


6064.225 881 


5 


7092.679 


7092.664 815 


6 


7614.622 


7614.603 506 


7 


8911.545 


8911.513 342 


8 


9095.696 


9095.679 497 


9 


10208.350 


10208.318 142 


10 


10869.289 


10869.255 077 


11 


11482.475 


11482.457 956 


12 


12353.799 


12353.766 422 


13 


12972.473 


12972.453 117 


14 


13690.455 


13690.436 602 


15 


14435.350 


14435.321 044 



TABLE XXII Eigenvalues in cm -1 of the Johnson asym- 
metric DWP consisting of the sum of a Morse and a 
Gaussian potentials V(r) — D[l — exp(—B(r — r a ))] 2 + 
Aexp(-C(r -r b ) 2 ). with A = 10 4 cm~ 1 ,B = 1.54A _1 ,C = 
200.0A" 2 ,D = 31250.0cm" 1 , r a = 1.5A,r b = 1.6A. John- 
son l|25r ) results are compared to the CFM. 



• Vibrational energy levels of Cold Molecules and ap- 
plication to the 23 Na2 molecule in the 0~ and l u 
electronic states. The Lennard- Jones molecule case 
is also studied. 

• Local and Non-Local Exchange problems. 

• Energy levels for regular and singular potentials 
(general ID case) 

Its mathematics is quite subtle since it enables one to 
find the Vibrational spectra of tenuous molecules where 
energies and distances are so remote from the ordinary 
short-range molecules case that special techniques should 
be developed in order to avoid numerical instabilities 
and uncertainties. 

The CFM has been tested succesfully in Long-range 
and Short-range potentials for Atomic and Molecular 
states and gives accurate results for bound and free 
states (f29l ;l3l"l;l35l). The tunable accuracy of our method 
allows to evaluate eigenvalues close to the ground 
state as well as close to highly excited states near the 
continuum limit to a large number of digits without any 
extrapolation. 

The CFM compares favorably with many different so- 
phisticated techniques such as those developed, for in- 
stance (to cite a few) , byRaptis et al. H^; |47[ ) , Johnson 
(IH), Hamilton et al. dS) and Rawitscher et al. dH). 



The CFM approach remains the same despite the wide 
variability of the mentioned problems. 

The VSCA integration method used gives the right 
number of all the levels and the variation of the eigen- 
value function F(E) definitely determines the total 
number of levels. Generally it requires performing 
analytically Taylor series expansion to any order of 
an arbitrary potential function that might require 
the combination of numerical, symbolic manipulation, 
functional fitting and analytic function continuation 
techniques (see refs (fiH li~7t)). Despite these challenges, 
the results of the VSCA are rewarding. 

The use of the RK4 and the VSCA methods jointly 
paves the way to a precise comparative evaluation of the 
accuracy of the spectra obtained. In practice, the RK4 
method can be used in optimisation problems whereas 
the VSCA is more adapted to the direct evaluation of 
the spectra. 

Since the CFM bypasses the calculation of the eigen- 
functions, it avoids losing accuracy associated with the 
numerical calculation specially with rapidly oscillating 
wave functions of highly excited states. This is specially 
needed in the study of the sensitive problem of determin- 
ing the vibrational spectra of cold molecules. 

We did not consider molecular rotation, nevertheless 
the CFM is adapted to solve accrately Ro- Vibrational 
problems as well as any RSE diagonalisation problem 
(see for instance ref. (|3ll)). 

Regarding spin dependent scattering problems we 
applied the CFM to local and non-local exchange 
integrating exactly the coupled equations originating 
from integro-differential equations. We benchmarked 
and tested our integration method at extremely small 
energies for momentum k > 10~ 8 and found that it 
performs again very well against stability error 65 and 
truncation error ex- Our results were validated with a 
highly accurate method, the recently developed S-IEM 
(|49|). The results we find (phaseshifts and scattering 
lengths) agree with the S-IEM results for the singlet, 
triplet cases for an arbitrarily small value of momentum 
k(k>l0~ s ). 

In conclusion, the (single channel and multichannel) 
CFM enables one to tackle precisely weakly bound 
states in atoms and (simple and cold) molecules as 
well as low-energy scattering of atoms and molecules 
and can be extended st raig htforwardly to Bose-Einstein 
condensation problems (|60h . 
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APPENDIX 
Atomic and other units 

In atomic and molecular physics, it is convenient to 
use the elementary charge e, as the unit of charge, and 
the electron mass m e as the unit of mass (although for 
some purposes the proton mass, m p , or the unified mass 
unit amu, is more convenient). Electrostatic forces and 
energies in atoms are proportional to e 2 /47reo , which has 
dimensions ML 3 T~ 2 , and another quantity that appears 
all over in quantum physics is h which has dimensions 
ML 2 T~ X ; so it is convenient to choose units of length 
and time such that 47re = 1 and % = 1. 

The atomic unit of length is then (by dimensional anal- 
ysis alone) 



m e (e 2 /47re ) 1 J 
I 



This is called the Bohr radius, or simply the bohr 
(0.529 A), because in the "Bohr model" the radius of 
the smallest orbit for an electron circling a fixed proton 
is (1 + m e /m p )aB- In the full quantum theory the par- 
ticles do not follow an orbit, but the expectation value 
of the electron-proton distance in the Hydrogen ground 
state is exactly (1 + m e /m p )aB)- 

The atomic unit of energy is the Hartree (27.2 eV) 
given by: 



47re h 



The unit of time is h/Eh- 

The Hartree is twice the ground state energy of the 
Hydrogen atom ^(l + mj/mp)" 1 ^ equal to the Ryd- 
berg (13.6 eV). In atomic and molecular spectroscopy, 
one uses rather the cm -1 an energy corresponding to a 
wavelength of 1cm or sometimes a frequency unit, the 
Hz. We refer the reader to the table below giving the 
conversion factors between the different energies. 





J 




eV 


Hz 


cm 1 


1 J 


1 




6. 24151. 10 18 


1.50919.10 33 


5.03411. 10 22 


1 eV 


= 1.60219.10" 


19 


1 


2.41797.10 14 


8.06547.10 3 


1 Hz 


= 6.62619.10" 


34 


4.13570.10" 15 


1 


3.33564.10 -11 


1 cm" 1 


= 1.96648.10" 


23 


1.23935.10 ~ 4 


2.99792. 10 10 


1 



Matching and boundary conditions 



and at r = oo: 



The CFM is based on the extraction of the eigenvalues 
from the zeroes of the eigenvalue function F(E) defined 
from the saturation of the left (r — > 0) and right (r — > 
+00) functions I -(E) and l+(E) given by the ratios of the 
canonical functions a(r) and /3(r). This was described 
previously in the case the boundary conditions are y(0) — 
2/(oo) = 0. In the general case we write: 

y(r) =y(r )a(E;r) + y'(r )/3(E;r) 
y'(r) = y(r )a>(E;r)+y'(r )[3>(E;r) (82) 

The canonical functions satisfy the conditions: 

a(E- r ) = 1, a'(E; r ) = 0, (3(E\ r ) = 0, 0'(E; r ) = 1 

(83) 

Let us rewrite the system [82] at the two boundaries 
r = : 



2/(0) 

y'(o) 



--y(r o )a(E;0) 
y(r )a'(E;Q) 



y'(r )l3(E;0) 
y'(r o )P'(E;0) 



(84) 



2/(oo) = y(r )a(E;oo) +y'(r )f3(E; 00) 
2/'(oo) = y(r )a'(E; 00) + y'(r Q )P'(E- 00) 

Extracting from above the left and right ratios: 

y'(r )' 



y(ro) 

y'(ro) 

Jj(ro) 



_ a(E;0)y'(0)-a'(E;0)y(0) 
- P'(E;0)y(0)-fJ(E;0)y'(0) 

a(E;oo)y f (oo) — a f (E-,oo)y(c>o) 
/3 / (£;oc)y(oc)- / 3(_E;oc)y'(oo) 



(85) 



(86) 



we can tackle several types of boundary conditions pro- 
viding below the energy functions I -(E) and l+(E) that 
determine the eigenvalue function: 



F(E)=l + (E)-l-(E) 



\y'(roY 






. y( r o) _ 


+ 


. 2/M . 



(87) 



Some of the boundary conditions are, for instance: 
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1. Standard (Dirichlct) case: 

The boundary conditions y(0) = y(oo) = yield: 



l-(E) = lim 



a(E;r) 



l+(E) = lim 



r^o /3(E;r) 
a(E;r) 



r^+oo (3(E;r) 



(88) 



2. Derivative asymmetric case type I: 

The boundary conditions y(0) = 2/(00) = yield: 



/_(£) = Km - 
l + (E) 



<*{E;r) 
""0 j3{E- r) 
a'(E;r) 



lim — - 

/?'(£; r) 



(89) 



3. Derivative asymmetric case type II: 

The boundary conditions y'(0) = 2/(00) = yield: 



l-(E) 
l + (E) = 



lim 



a'(E;r) 
■""a ' P'(E;r) 

lim -^M 

(3(E; r) 



(90) 



4. Derivative (Neumann type) boundary conditions: 
The boundary conditions j/(0) = y'(oo) = yield: 



I _(E) = lim - 
l + (E) = 



a'(E; r) 
""0 /3'(E;r) 
a'(E;r) 



lim — „ 

(3'(E;r) 



(91) 



These formulas can also be generalised to arbitrary 
(Cauchy or mixed type) boundary conditions: 



a lV (0) + b lV (Q) =ci 
a 2 y(oo) + b 2 y(oo) = c 2 



(92) 



and to the multichannel case. For instance, in the ID 
symmetric Double Gaussian potential the boundary con- 
ditions correspond to above case No. 3 whereas the asym- 
metric case corresponds to case No. 2. 
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